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NUMERICAL PROCEDURES FOR CALCULATING REAL FLUID 
PROPERTIES OF NORMAL AND PARAHYDROGEN 
by Fredric N. Goldberg and Angela M. Haferd 
Lewis Research Center 

SUMMARY 

Numerical procedures are described for calculating real fluid properties of normal 
or parahydrogen. The library of single function calls can be used efficiently without ini- 
tial estimates. When physical conditions are known, engineering estimates of density 
may be included for additional speed in calculation. 

Each of the procedures is designed to calculate a specified function of one or two in- 
put parameters. Independent variables consist of the combinations temperature- density, 
pressure-density, temperature-pressure, enthalpy-pressure, or entropy-pressure. 

Ideal and saturation properties may also be calculated as functions of either temperature 
or pressure. 

The program is written in separate modules that can be modified easily in the future 
to include new data. These modules can be used as a subset of functions when computer 
space is not available for the entire library in complex programs. 

The discussion includes techniques used, resolution of inherent physical and mathe- 
matical problems, and consistency of results based on National Bureau of Standards 
data. Computer time and storage estimates are given with a manual of instructions for 
the programmer. A listing of the program is included as used in the library of functions 
at the Lewis Research Center. 


INTRODUCTION 

Present research effort in advanced propulsion systems involves the use of molecu- 
lar hydrogen as both a propellant and a coolant. The performance analysis of such sys- 
tems as regeneratively cooled rocket engines requires accurate knowledge of transport 
and thermodynamic properties in both the liquid and the vapor phase. 

Several methods have been published (refs. 1 to 4) for specific applications such as 



theoretical analysis (ref. 3) or table computation (ref. 4). In the procedures used in ref- 
erence 3, close initial estimates are important for convergence. These predictions can 
be difficult for a user to furnish, and an incorrect estimate can cause the calculation to 
go out of bounds or iterate indefinitely. In tabular methods (ref. 4), inconsistencies re- 
sult from interpolations across phase boundaries. Multiple function calls used in refer- 
ence 4 for increased speed are not in a form suitable for writing complex engineering ex- 
pressions where single functions are designated as parameters. 

In reference 2 the correlations proposed for recent National Bureau of Standards 
(NBS) tables are extremely difficult to program and require more computer storage than 
is available when library systems are used. The equation of state proposed in refer- 
ence 1 by NBS provides a high degree of consistency except in the critical region and for 
temperatures above 400° K. With modifications and additions in these areas, this equa- 
tion of state can be used effectively. Related correlations from references 1, 2, 5, and 
6 were incorporated in the procedures described in the present report. 

The library of functions presented by the authors is currently in general use at the 
Lewis Research Center for normal or parahydrogen, as required in actual use. The tem- 
perature ranges from 14° to 1500° K, and the pressure ranges from 0. 1 to 340 atmos- 
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pheres (approximately 0. 1x10 to 340x10 N/m ). Typical results for all functions ex- 
cept specific heat give a consistency of the order of 0. 02 to 0. 05 percent in comparison 
with tabulated NBS data. 

The method used herein follows systems programming practice to simplify effort for 
the user. A single function is specified, and no initial estimate is required. Phase 
boundaries are automatically defined in all the procedures. The library is in the form of 
separate modules that can be updated as new data are obtained. 

This report presents all pertinent material used in developing the programs of the 
library. Instructions for use of the functions are given in appendix A. Although the user 
merely needs to specify the functions as listed in actual application, the report is written 
to provide the scientific programmer with equations and derivations for reference pur- 
poses. For those interested in programming techniques and the solution of the physical 
problems, these details of the system are developed. The FORTRAN IV program listing 
in appendix B gives further details for a closer examination and provides coefficients and 
constants that are not in the text. 


A, B, C, D, F constants 



SYMBOLS 


2 


sonic velocity 

specific heat at constant pressure 



C g £ specific heat of saturated liquid 

C y specific heat at constant volume 

H enthalpy 

K thermal conductivity coefficient 

m molecular weight 

P pressure 

P crit critical pressure, 12.77 atm; 1. 29392X10 6 N/m 2 

R gas constant 

S entropy 

T temperature 

T crit critical temperature, 32. 984° K 

V volume 

Z compressibility factor, P/pRT 

y specific heat ratio, Cp/C y 

p shear coefficient viscosity 

p density 

P cr it critical density, 0. 0152672 (g) (mole) /cm 3 ; 30. 7744 kg/m 
a defined by eq. (3) 

Subscripts: 
crit critical value 

cs program constant 

h high 

id ideal gas 

l low 

lim limit 

liq liquid 

o base line reference 

sat saturated line 

si saturated liquid 



sv saturated vapor 
x mixture under vapor dome 


PRESSURE-TEMPERATURE-DENSITY PROCEDURES 
Basic Equations 


Equations of state are used to represent the basic physical relation among density, 
pressure, and temperature. 

A modified Benedict-Webb-Rubin equation of state (ref. 1) is used for temperatures 
up to and including 400° K: 


o / a 4 a e \ q 

P = a 1 pT + p [ a x a 2 T + a * + -± + — \ ) + P 


T t 2 T 4y 
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The coefficients to a 1? are specified in two groups for regions designated as 
either vapor or liquid. Combinations of temperature and density identify the two regions. 

For temperatures above 400° K, a high-temperature equation of state from refer- 
ence 5 is used: 


P = RTp exp 


B(T)p + C(T)p 2 ] 


( 2 ) 


where 


B(T) = BjT"l/ 4 + B 2 T" 3 / 4 + BgT' 5 / 4 


and 


C(T) = C 


t-3/2 

l 1 


+ C 2 T " 2 


Equations derived from reference 1 are used to evaluate the two-phase region. 



Equations (3) and (4) define the liquid boundary: 


p sat = p crit - t^)l 3 E b n[' r(T) ! n 

n=0 

4 

<**) - E a n< T orit - T > n/3 

n=l 


( 3 ) 


Psi = Pent + a < T > ( 4 > 

The saturated vapor boundary is defined as the values of density at the intersection of 
isotherms (eq. (1)) with saturation pressure (eq. (3)). An ideal gas is defined as one that 
is described by the perfect gas law, where P = RpT, but has a temperature -dependent 
specific heat. 


Region Determination and Initialization 

All property evaluations involving equations of state must be preceded by a region- 
finding process to determine which equation of state applies. The region-finding process 
also provides additional information required to initialize iterations and to specify bound- 
ary values. These regions are described in the diagrams of figure 1. 

Since all state equations used give pressure as a function of temperature and density, 
the combination of temperature and density as independent variables directly determines 
the region. For the cases where independent variables are pressure and density, or 
pressure and temperature, indirect methods of region determination have to be made. 
These implicit determinations depend on the property that both derivatives (3P/ 3T)^ and 
(3P/3p)rp are nonnegative. 


Calculation Procedure 

The procedures for evaluating the pressure-density-temperature surface are pro- 
grammed in the form 


P(p, T) pressure 
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T(P,p) 


temperature 


p (T, P) density 


The evaluation of pressure is a direct substitution of the values of density and tempera- 
ture in the appropriate state equation. 

Temperature and density procedures involve a method of successive approximations 
in which a sequence of trial values is generated by a Newton type iteration: 


n+1 


= T 


n 


P(T 

3P 

dT 


n’P> - P 
(T n >P) 


p n+l p n 


P(T,P n ) - P 

f (T’Pn) 

dp 


Indicated partial derivatives are formal evaluations of derivatives of the state equation 
applicable for the region. 

The values of iterants are governed by algorithms that recognize the behavior of the 
particular function and its derivatives in each region of interest. This technique prevents 
values from exceeding the range indicated by the region finding process and guarantees 
convergence. 


Critical Region Calculations 

In the vicinity of the critical region, the density calculation is extremely unstable. 
Large errors occurred in calculating functions and derivatives of the equation of state. 

To eliminate this problem, an empirical correlation was substituted. Two sets of density 
values p^ and p^ were fit empirically for isotherms in the temperature range between 
32. 5° and 35° K. The variation of density with pressure in this region was then de- 
scribed by an interpolation formula. The criteria used to isolate the region are symbol- 
ically shown in figure 2. 

The resulting representation used for the density procedure is summarized as fol- 
lows: 
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For 32. 5 < T < T 


crit' 


PGy T> < p < P sat< T > 

^sv " ~ P sat^ + 

' [ P saf P ^ T )] SV 

P sat< T > < p * P H’ T > 

( Ph-' 3 slX p - p <Ph> T )1 ( 

- P satl +Ph 


For 


T crit — T < 35 L 


K: 


P(p z ,T) <P<P(p h ,T) 


(P h “ P?)[ p " p (Ph> T )] 

P = — — + p h 

[P0p h ,T) - Pfy,T)] 

State equation (1) is used to calculate the values P (p^, T) and P(p h , T) with the use of 
liquid or vapor coefficients as indicated. The result of this technique is a distortion of 
the vapor-liquid dome to avoid numerical difficulties and errors. 


THERMODYNAMIC AND TRANSPORT PROCEDURES 
Functions of Temperature and Pressure 

The first group of four procedures for thermodynamic and transport properties is 
programmed in the form generally used in data reduction. Pressure and temperature are 
the given parameters: 

H(P, T) enthalpy 

S(P, T) entropy 

C (P, T) specific heat, also includes C y and C Q 
p(P, T) viscosity, also includes thermal conductivity K 
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In each of these procedures, density p(P, T) is first calculated, which automatically 
identifies the region. Then the property is calculated by combining a function of temper- 
ature only with an isothermal deviation term. (See appendix C for details. ) The devia- 
tion term reflects the real gas effects, and it is represented by either functions of the ap- 
propriate state equation or other empirical correlations. These relations are summa- 
rized as follows: 



C 


o 




1/2 


C 


v 


M = M q (T) + p 1 (T,p) 

K = K o (T) + K t (p) 

All derivatives and integrals in these equations are evaluated in closed form from the 
equations of state. 
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Functions of Pressure and Enthalpy or Entropy 

Two procedures for calculating temperature and density for a given value of pressure 
with enthalpy or entropy are of the form T(P, H) and T(P, S). 

Since the equations of state are represented as explicit functions of p,T, it is pos- 
sible to express H or S as a function of the same two parameters. These calculations 
involve simultaneous solution of two nonlinear equations. To guarantee a convergent se- 
quence, the problem was linearized. A sequence of temperatures is generated with the 
use of a one-dimensional Newton type iterative scheme. An intermediate calculation of 
density is performed for each temperature and pressure. Initial values of temperature 
are chosen through an algorithm that avoids crossing region boundaries. Infinite loops in 
the computation are thus avoided for the intermediate iterates. 

The initial region calculation in this case may determine that the state point (H, P) or 
(S, P) lies within the two-phase region. When this occurs, the problem is reduced to cal- 
culating the saturation temperature corresponding to the given pressure. The density for 
the mixture region is then calculated from the following equations that represent an aver- 
age density variation with mixture quality: 


p _ p sv^sZ 

x [xp si + (1 - *>P SV ] 


where 


H - H 


x = 


Sl 


H sv - H si 


or 


x = 


S - 


sv 



Functions of One Parameter 

To facilitate the programming of data reduction problems in which only one param- 
eter can be measured and a given phase value of the fluid is to be assumed (such as sat- 
urated liquid, saturated vapor, mixture quality, etc. , ), procedures are included to cal- 
culate saturated or ideal values of the following properties as functions of temperature: 
P ga t(T) is given by equation (3), p g ^(T) by equation (4), P SV (T) occurs at the intersection 
of equations (1) and (3), 
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C P, id< T > = A id + B id T - c id T + D id T + E id T + F id T ' 
H id (T) - H id (T 0 ) + £ C p; id (T)dT 


*uF> = S ld< T «> 


c/rp 


5 xt / 2 \ 


where S id (T Q ) is evaluated at a pressure of 1 atmosphere (1. 01x10 N/m ) 


p . /» p sv r 

H SV (T) = H ,(T) + 55* _ + / J 

sv ia p - RT / 

Psv ./) P 


P T /9P1 


2 2 \0T 

P P x 'P 


dp 


S sv (T) = S id( T) " R ln ( RT P s J + 


SV 7 


/ 


sv 


R _ J_ /3P\ 

P r 2 WA 
F P X / D 


dp 


V T > = W - / T c «< T > dT - / 


sZ 


dP 


(T ) P 
sZ u o ; 


where 


V T > ■= 


A T 
cs 


< T crit ~ T ) 


1/2 


+ B + C T + D T 2 + E T 3 
cs cs cs cs 


(5a) 


(5b) 


S sz< T) = W + 


C sZ ( T ) 



dT 


C v (T) = C 
sv 


P, id 


(T) - R + 




dp 
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\ (T) = C(p 1 ,T) + 
si 


/ 


p sZ /9C 



where C (p 1 ,T) is heat capacity at constant density, as tabulated in reference 2, and 

V X o 

Pj = 0.037821 gram-mole per cubic centimeter (or 76.237 kg/m ). The quantities Cp 


or C v 


and C are calculated with the use of relations summarized for functions of 


si 


sv 


temperature and pressure. However, caution must be exercised in the use of Cp and 
C D at the saturation line, because the values returned are calculated as limits as the 
density approaches saturation but do not reflect the physical fact that the actual value of 
the property may become unbounded at the onset of cavitation or two-phase flow. 

In addition to the given functions of temperature, a procedure is included for calcula- 
tion of saturation temperature based on an input pressure. 

The quantity T ga ^.(P) is obtained by an iterative calculation based on an empirical 
equation from reference 3: 


P sat = ex P ^sat) 


(6) 


where 


F <T Sa t> - A ts + - Bt f - c ~ + D ts T sat 
1 sat + Ss 

EFFECT OF COMPOSITION ON PROPERTIES OF NORMAL OR PARAHYDROGEN 

The primary use of the library of property functions is for parahydrogen composition. 
An option is included for normal hydrogen composition in the vapor region in the proce- 
dures for enthalpy, entropy, specific heat, and transport properties. When the composi- 
tion is 75-percent orthohydrogen and 2 5- percent parahydrogen, the gas is considered to 
be in the normal mode. In the liquid region, the presence of normal composition would 
be extremely rare in actual application and is not considered. 

The pressure-density-temperature surface for hydrogen is nearly independent of 
orthopara percentage. State equations, real gas deviation, and excess function can be 
treated as independent of the effect of composition. However, the major difference in 
property values is in the base line or ideal gas value. In the procedures for enthalpy, 
entropy, and specific heat, provision is made to calculate the base effects for either nor- 
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mal or paracomposition as designated by the user. 

The effect of composition on transport properties is significant only in thermal con- 
ductivity. The proper value of the ideal gas specific heat is used to reflect the difference 
in composition between normal and para. In the viscosity function, the effect is not sig- 
nificant but the programming of the two functions is related in the procedure call, and the 
option is used for normal or para. 

Liquid hydrogen is considered to be equivalent to para in the research applications in 
heat transfer for which the functions are used. However, a relatively good estimate of 
liquid-hydrogen normal functions can be derived by correcting values for parahydrogen 
by the difference in ideal gas values. The library of functions includes ideal gas proper- 
ties for this purpose and for other uses. 


USE OF PROCEDURES AND RESULTS 
Consistency and Validity of Results 

Error analysis is considered in two categories. The first is internal consistency 
within the group of functions, and the second is comparison of results with reference 
values or source data. Physical and mathematical consistency is considered most im- 
portant in the experimental and analytical computations. Misleading physical results 
must be prevented, such as an isobaric increase of enthalpy indicating a temperature 
drop or density rise. The mathematical consistency among the related functions and 
their inverses is kept to within approximately 0. 1 percent. Slightly larger errors occur 
near the critical pressure-temperature point that must be recognized in any application. 
Certain tolerances are also used at the saturated liquid boundaries for mathematical rea- 
sons. 

Comparison of results obtained for pressure-density-temperature relations with ref- 
erences 2, 5, 7, and 8 shows agreement to within approximately 0. 02 percent for all re- 
gions except critical. The error of approximation compares with the experimental ac- 
curacy represented in the available sources of data. In some functions, such as enthalpy 
and entropy, the agreement within the published sources varies slightly because of refer- 
ence conditions. The specific-heat function is uncertain in the saturated region, and near 
the critical point it becomes unbounded. This function changes rapidly at other low pres- 
sures, as shown in figure 3. However, overall results compare well at all points, except 
the critical region, with errors within 0. 5 to 2. 0 percent. 

The agreement for transport properties is completely consistent with recent informa- 
tion (refs. 3, 6, 9, and 10) and the same mathematical method is used with identical em- 
pirical relations. However, the validity of referenced correlations for viscosity and 
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thermal conductivity is not definitely established at present. Computer-generated plots 
of the properties of parahydrogen as represented in this report are included in figures 3 
to 8. 


Execution Time for Function Calls 

For programming applications, the average time estimate for each function call in- 
volving state equations and iteration is at the rate of approximately 200 per second. Cal- 
culation speed can be increased by a factor of 3 to 4 when close estimates of density are 
furnished in any function of temperature and pressure. Successive calculations of func- 
tions in which there are only small variations in density also automatically provide in- 
creased speed. When density is calculated along with another function, it is also pro- 
vided as a secondary result in common. Other related functions are also grouped that 
can be computed simultaneously, such as viscosity- conductivity, or a combination of spe- 
cific heats Cp and C v , and sonic velocity C Q . 

Explicit functions, such as pressure in state equations, or saturated region relations 
and ideal functions, are calculated at a rate of several thousand per second. 


Programming Use of Library of Functions 

The list of function calls with instructions for programming is given in appendix A. 
The prefix identifying the library represents the system of units being used. The letters 
BW refer to the internal library units corresponding to the Benedict-Webb-Rubin equation 
of state (eq. (1)), which is in the calorie-gram-centimeter system. The particular func- 
tion is specified by a number according to the listed parameters given. 

For input and output in British thermal units, the prefix BT is used in the function 
call. The prefix BI is used for input-output in the SI system. The function number is 
also modified to indicate that the composition of the hydrogen is normal rather than para- 
hydrogen. 


CONCLUDING REMARKS 

The procedures presented have been applied at the Lewis Research Center effectively 
in various experimental and theoretical research projects. The applications included 
large programs, such as convective cooling analysis, reduction of data from nuclear en- 
gine simulations, rocket engine design and performance calculations, and correlations for 
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experiments with pressurized tank storage and pumping of hydrogen. 

Evaluations of performance parameters have been made efficiently and smoothly 
proving the consistency and versatility of the functions used. In comparison with National 
Bureau of Standards tables, the properties calculated were more consistent than those ob- 
tained with the use of other methods. Improvements were made especially in the critical 
region and in the liquid- vapor phase boundary definition. 

The techniques originated by the authors provide solutions for similar problems that 
would be encountered in applying state equations for properties of other fluids. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, September 13, 1967, 

125-23-02-10-22. 
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APPENDIX A 


PROGRAMMING INSTRUCTIONS FOR HYDROGEN PROPERTIES 

Calling Sequence 

The programming function call is written in the form of either BW(N, PI, P2), 

BT(N, PI, P2), or BI(N, PI, P2). The prefix designates which system of units is to be 
used in the conversion: BW is used for cgs units (calorie- gram- centimeter- second) sys- 
tem, BT is used for British thermal units system, and BI is used for the SI system. 

The integer N specifies which thermodynamic property is to be computed. Param- 
eters PI and P2 are input combinations depending on N, as shown in table I. 

Related output is also given in NAMED COMMON from certain procedures. 


TABLE L - VALUES OF N a FOR SELECTION OF THERMODYNAMIC PROPERTY 
[For British thermal units, use prefix BT in place of BW; for SI units, use prefix BI.] 


Function 

call 

Designated function 

Related output 

BWX 

BWD 

BWCVGL 

BWCOND 

BWCO 

BW(1, T, P) 

Density, D 

X 

D 

— 


— 

BW{2, T, P) 

Enthalpy, H 

X 

D 



— 

BW(3, T, P) 

Entropy, S 

X 

D 

--- 


— 

BW(4, H, P) 

Temperature, T 

X 

D 

--- 


— 

BW(5, S, P) 

Temperature, T 

X 

D 



— 

BW(6, T, P) 

Specific heat, CP 

X 

D 

cv 


CO 

BW (7, T, P) 

Viscosity, \±, xlO 4 

X 

D 

— 

KxlO 4 

— 

BW(8, T, D) 

Pressure, P 

X 

— 

— 


— 

BW(9, D, P) 

Temperature, T 

X 

— 

--- 


— 

BW(10,T, O) 

Saturation pressure, P 

--- 


— 


— 

BW(11, P, O) 

Saturation temperature, T 

--- 

--- 

--- 


— 

BW(12,T,0) 

Density of saturated liquid 

— 

--- 

— - 


— 

BW(13, T, O) 

Density of saturated vapor 



--- 


— 

BW(14,T,0) 

Specific heat of ideal gas 


— 

--- 


— 

BW(15, T, O) 

Specific heat of saturated liquid 

— 

D 

cv 


CO 

BW(16,T, O) 

Specific heat of saturated vapor 


D 

cv 


CO 

BW(17,T,0) 

Enthalpy of ideal gas 

--- 


--- 


— 

BW(18, T, O) 

Enthalpy of saturated liquid 

--- 


-- - 


--- 

BW(19, T, O) 

Enthalpy of saturated vapor 


D 

— 


— - 

BW(20, T, O) 

Entropy of ideal gas 

— 

— 



— 

BW(21,T, O) 

Entropy of saturated liquid 

— 

D 

--- 


--- 

BW(22,T, O) 

Entropy of saturated vapor 


D 

— 


— 


a For normal hydrogen, N is replaced by N + 200. 
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Physical Units 


In table n the physical unts and symbols used in the various function calculations are 
presented. 


TABLE n. - SYSTEM OF UNITS FOR VARIOUS PROPERTIES 


Property 

Symbol 


Physical units 




BW 

BT 

BI 

Temperature 

T 

°K 

°R 

°K 

Pressure 

P 

atm 

lb/sq in. 

N/m 2 

Density 

D 

g/cc 

lb/cu ft 

kg/cu m 

Enthalpy 

H 

cal/g 

Btu/lb 

J/kg 

Entropy 

S 

cal/(g)(°K) 

Btu/(Ib)(°R) 

J/(kg)(°K) 

Specific heat at con- 

c p 

cal/(g)(°K) 

Btu/(lb)(°R) 

J/(kg)(°K) 

stant pressure 




Specific heat at con- 

c v 

cal/(g)(°K) 

Btu/(Ib)(°R) 

J/(kg)(°K) 

stant volume 





Sonic velocity 

c o 

cm/sec 

ft/sec 

m/sec 

Viscosity 

g/(cm)(sec) or PxlO^ 

[ Ib/(ft) (sec) ]x 10 4 

[(N)(sec)/sq m]xl0^ 

Thermal conductivity 

K 

[cal/(cm)(sec) (°K)]xl0^ 

[Btu/ (ft) (sec) (°R)]xl0 4 

[J/(m)(sec)(°K)]xl0 4 

Number designating 
region 

X 

Liquid (-1), vapor (2), or 
liquid- vapor mixture (0-1) 




In the liquid- vapor two-phase region, x is equal to the quality of the fluid. 


Applicable Range of Validity 


The range for which the function calls are valid can be summarized in terms of tem- 
perature and pressure approximately as follows: 


Temperature: 

°K 

°R 

Pressure: 
atm .... 
psi . . . . 
N/m 2 . . . 


. ... 14 to 1500 
. ... 25 to 2700 

. . . . 0. 1 to 340 
. . . . 1 to 5000 

10 000 to 340X10 5 
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Intermediate Results in NAMED COMMON 


Depending on the choice of N, the following properties result in common: 

BTX, BWX, or BIX quality or region 

BTD, BWD, or BID density 

BTCVGL, BWCVGL, BICVGL specific heat at constant volume 

BTCOND, BWCOND, BICOND thermal conductivity 

BTCO, BWCO, or BICO sonic velocity 

NAMED COMMON is required in the main program whenever any of this related out- 
put is to be used, and is declared in double precision as follows: 

COMMON/BWDXCK/BWD, BWX, BWCVGL, BWCOND 

COMMON/BTDXCK/BTD, BTX, BTCVGL, BTCOND, BTCO 

COMMON/BWSONV/BWCO 

COMMON/BIDXCK/BID, BIX, BICVGL, BICOND, BICO 
DOUBLE PRECISION BTD, BTX, BTCVGL, BTCOND, BTCO 
DOUBLE PRECISION BWD, BWX, BWCVGL, BWCOND, BWCO 
DOUBLE PRECISION BID, BIX, BICVGL, BICOND, BICO 

Normal Hydrogen 

For properties of normal hydrogen, the value of N given in table I is replaced by 
N + 200 in all cases where the region is applicable as vapor. If the programmer speci- 
fies normal hydrogen in the function call, the value returned when the region is not found 
to be vapor will be for parahydrogen composition. 

Saturation Conditions 

If saturation conditions are determined to exist where pressure and temperature are 
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the independent variables, the results returned are those corresponding to saturated liq- 
uid conditions. 


Special Properties 

In any function that uses temperature and pressure as independent variables, an en- 
gineering estimate of density may be provided. This value is put in BWD, BTD, or BID 
depending on the system of units used. 

In general, a single result is returned as the value of an expression. For example, 
in computing the free energy H - TS given P and T, the programmer writes 

F = BW(2, T, P) - T X BW(3, T, P) 

In this case, the intermediate results H = BW(2,T, P) and S = BW(3,T, P) are not auto- 
matically stored. 


Computer Storage Requirements 

Approximately 9000 (decimal) computer stores are required. 
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APPENDIX B 


PROGRAM LISTINGS 
Grouping of Related Modules 

Group (1) BW program used for cgs function calls 

Group (2) Basic P, p, T group (requires module of group 1): 

BWDENS, BWBLOK, LIQUID, VAPOR, 

PK, PPK, COEF, BWPSAT, BWROSV, 

ROKSL, BWTSAT, DPDT, BWTEMP, BWPRES 

Group (3) Thermodynamic group of H, S, C p functions (requires modules of groups 1 
and 2) : 

PARALO, PARAHI, NORMLO, NORMHI, COE FID, 

BWCPIG, BWHIDG, BWSIDG, BWSSTL, BWHSTL, 

BWSSV, SKSL, BW SSL, HKSV, HKSL, SKSV, 

HDEV, BWHSV, BWHSL, SDEV, BWENHY, 

BWERPY, BWTOFH, BWTOFS, BWFD2, BWCPGL 

Group (4) Transport properties p, K for all temperatures and high-temperature func- 
tions P, p, T, H, S, Cp, C y , Cq (requires modules of groups 1, 2, and 3): 

BWVTSC, HTCPIG, HTHIDG, HTSIDG, HTHDEV, 

HTENHY, HTSDEV, HTERPY, HTTOFS, HTDPDT, 

HTDENS, HTFD2, HTTOFH, HTCPGL, HTTMP, HTPKS 

Group (5) BT program used for British thermal units function calls (optional) 

Group (6) BI program used for SI units function calls (optional) 

Group (7) Function NT AG (optional for omitting tagged data check) 
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FORTRAN IV Program Listing For Library of Functions 


IBFTC BW DECK 

SUBPROGRAM USED FOR ALL PROPERTIES IN BW UNITS (CGS) SYSTEM 
FOR PARAHYDROGEN N=1 - 22 OR NORMAL HYDROGEN N=2Cl-222 
IN CALOR I ES-GRAM- SECOND SYSTEM. P=.l - 340 ATM 
TEMPERATURES 14 TO 1500 DEG KELVIN 
FUNCTION 8W(N,Pl t P2) 

COMMON/ F I XED/MODEt ITERS 

COMMON/ BWDXCK/ BWD» BWX t BWCVGL * BWCOND 

COMMON/BWSONV/BWCO 

DOUBLE PRECISION BWD, BWX , BWC VGL , BWCOND 
DOUBLE PRECISION BWCO 
MODE = N/100 
T 400=400 • 

M= N— 100#M0DE 
BW= P1 + P2 

IF (NTAG(BW) .NE.O) GO TO 333 

GOTO (l t 2 r 3,4, 5, 6,7, 8, 9,10*11*12,13, 14 ,15,16, 17, 18, 19, 20, 21, 22 ),M 

333 BWD= BW 
BWX=BW 

IF (M.NE.6) GO TO 334 

R WCVGL=B W 

BWCO-BW 

RETURN 

334 IF (M.NE.7I GO TO 335 
BWC0ND=6W 

335 RETURN 

1 I F ( P 1 .GT . T 400 ) GO TO 101 
B W= BWDENSf P1,P2> 

RETURN 

101 BW= HTDENS(P1,P2) 

RETURN 

2 IF( P1.GT.T400) GO TO 102 
BW= BWENHY(P1,P2) 

RETURN 

102 B W = HTENHY(P1,P2> 

RETURN 

3 IFfPl.GT.T 400 ) GO TO 103 
BW- BWERPY(P1,P2) 

RF TURN 

103 B W= HTERPY ( P 1 , P2 ) 

RETURN 

4 H= BWENHY ( 400 . , P2 ) 

IF(Pl.GT.H) GO TO 104 
BW= BWT0FH(P1,P2) 

RETURN 

104 BW= HTT0FH(P1,P2) 

RETURN 

5 S- BWERPY ( 400. , P2 ) 

IF(Pl.GT.S) GO TO 105 
BW= BWTOFS ( PI , P2 ) 

RETURN 

105 BW= HTT0FS(P1,P2) 

RETURN 

6 IFfPl.GT.T 400 ) GO TO 106 
BW= BWCPGL(P1,P2) 

RETURN 

106 BW= HTCPGL l P 1 , P2 ) 

RETURN 

7 BW= BWVISC(P1,P2) 

RETURN 
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8 IF(P1.GT.T 400 ) GO TO 108 
BW=BWPRES(P1,P2) 

RETURN 

108 BW-HTPRES ( Pit P2 ) 

RETURN 

9 P = BWPRES(T 40 0 , P 1 ) 

I F { P2 .GT • P ) GO TO 109 
BW=BWTEMP(P1,P2) 

RETURN 

109 BW=HTTEMP(P1,P2) 

RETURN 

10 BW= BWPSAT(Pl) 

RETURN 

11 BW= B WT SAT { P 1 ) 

RETURN 

12 BW= ROKSL (PI ) #2 • 0 1 57 
RETURN 

13 BW= BWROSV(Pl) 

RETURN 

14 I F { P1.GT.T400) GO TO 114 
BW= BWC P I G ( P 1 ) 

RETURN 

114 BW= HTCPIG ( PI ) 

RETURN 

15 P S A T = BWPSAT (PI ) + . 0 0 1 
BW= BWCPGL(PltPSAT) 

RETURN 

16 PS AT= BWPSAT(P1)-.001 
BW= BWCPGL ( PI ,PSAT ) 

RETURN 

17 I F ( P 1 . GT • T 400 ) GO TO 117 
BW= BWHIDG(PI) 

RETURN 

117 BW» HTHIDG(PI) 

RETURN 

18 BW= BWHSTL (PI ) 

RETURN 

19 PS AT = BWPSAT (PI ) + . 001 
B W = BWHSV(PSAT) 

RETURN 

20 I F ( P 1 • G T • T400 ) GO TO 120 
BW= BUS I JG(P1) 

RETURN 

120 BW= HTSIDG(Pl) 

RETURN 

21 BW= BWSSTL(Pl) 

RETURN 

22 PS AT = BWPSAT(Pl) 

BW= BwSSV ( PSAT ) 

RETURN 

END 

IRFTC BWDENS DECK 

CALCULATES DENSITY IN C.G.S. UNITS AS FUNCTION 
OF TEMPERATURE AND PRESSURE. REQUIRES BASIC MODULES BWBLOK, 
BWPSAT , VAPOR, LIQUID, PK, PPK, COEF, AND ROKSL 
DOUBLE PRECISION FUNCTION BWDENS (T6,P6) 

COMMON/ FI X ED/ MODE , I TERS 
COMMON/ B WDXCK/ BWD,BWX, BWC VGL,BWCOND 
COMMON/ PTD/ P,T,D,ROK, DELTA, XV AL 
COMMON/ WORK S / WK S , W 1 , W2 , W3 , W4 , S 



COMMON /KUNS/TCRIT,PCRIT f ROCR IT t TOL » DAI R 
DOUBLE PRECISION T6,P6,BWPSAT,R0KSLtPK,PPK 
DOUBLE PRECISION BWD,BWX ,BWCVGL ,BWCOND 
DOUBLE PRECISION P, T , D, ROK, DELTA, XVAL 
DOUBLE PRECISION WKS , W1 , W2 ,W3 ♦ WA ,S 
DOUBLE PRECISION TCR I T , PCR I T * ROCR I T » TOL , DAI R 
LOGICAL ESTD 
ESTD= * TRUE • 

ROK=BWD/2 .0157 
R0ID=P6/T6/82.082 
IF (TCRIT-T6) 2*1,1 

1 WKS = RWPSAT ( T 6 ) 

CALL LIQUID 

IF (WKS - P6 ) 3,3,22 
3 IF ( P6 - WKS - .0001) 5, A, A 
A XVAL = -1 . 

B0UNDL=R0CR IT + S 
UB0UND=.05 
RESTR= .05 

IF (T6.GT.32.5) GO TO AA 
GOTO 10 

5 ROK = ROKSL (T6) 

CALL COEF ( T6 ) 

WKS= PK(ROK) 

XVAL » 0 

GOTO 12 

2 CALL VAPOR 
BUUN DL = 0 • 

UBOUND= .0 5 
XVAL a 2 

IF (T6.LT.35.) GOTO 66 
IF ( P6 - PCR I T ) 6,6,7 

6 RESTR=ROID 
GO TO 10 

7 RESTR=.05 

10 CALL COEF (T6) 

15 DO 11 I T ERS= 1,10 

IF (ESTD) GO TO 1000 

16 X = PK(ROK) 

WKS = P PK ( ROK ) 

D E L T A = ( X - P6J/WKS 

ROK=ROK - DELTA 

I F ( DABS ( DELT A ) - TOL) 12,12,11 

11 CONTINUE 
GO TO 12 

1000 IF ( ROK.GT.BOUNDL.AND.ROK.LT .UBOUNDIGO TO 16 
ROK = RESTR 
ESTD= .FALSE. 

GO TO 15 
60 UBOUND=ROLOW 
RESTR=ROI D 
GO TO -15 
22 CALL VAPOR 
8UUNDL=R0ID 
UBOUND=ROCRIT 
RESTR=ROI D 
ROK=RO I D 
XVAL- 2. 

ROLOW= — . 0A3 + .0016* T6 
ROH I GH= ROCRIT - S 
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PHI GH= WKS 
CALL COEF ( T6 ) 

P LOW = PK(ROLOW) 

I F ( T6.LT.32.5) GO TO 15 
UBOUND=ROLOW 
IF { P6 . L T • PLOW ) GO TO 6 
GO TO 100 

44 ROH I GH= .086 - .002 * T6 
RO'JNDL=ROH I GH 
ROL 0 W= S + ROCRIT 
PLOW = WKS 
CALL COEF { T6 ) 

PHI GH= PK(ROHIGH) 

IF (P6.GT.PHIGH) GO TO 10 
GO TO 100 

66 ROLOW= — .043 + .0016* T6 
ROH I GH = .086 -.002 *T6 
CALL COEF ( T6 ) 

PLOW= PK(ROLOW) 

IF ( P6.LT.PL0W) GO TO 60 
CALL LIQUID 
CALL COEF ( T6 ) 

PHI GH = PK(ROHIGH) 

BOUNDL= ROHIGH 
RESTR= .05 

IF ( P6.GT.PHIGH) GO TO 15 

100 ROK = ROI!IGH-{ PHIGH-P6 ) / ( PHIGH-PLOW )* ( ROH IGH-ROLOW ) 

P H I G H = PK(ROK) 

12 0= ROK * 2.0157 
BWDENS = D 
BWD = D 
R W X= XVAL 
RETURN 
END 

IBFTC BWBLOK DECK 

PROGRAM CONSTANTS FOR CRITICAL TEMPERATURE, PRESSURE, DENS I TY, 
COEFFICIENTS FOR EQUATION OF STATE (1) 

BLOCK DATA 

COMMON/ FIXED/MODE, ITERS 

COMMON/ BWDXCK/BWD,BWX,BWCVGL,BWCOND 

COMMON/ PTD/P,T,D, ROK, DELTA, XVAL 

COMMON/ PVS/PP,BT,CT,BP,CP,H,D1,D2,D3,D4,EXPF 

C OMMON/CVS / CS,CH,F,E,U,C,B,AA 

COMMON/ FTS/AF,BF,CF,DF,EF,FF,B1,C1,C2,CK,E1,EK 
COMMON/ WORKS/ WKS ,W1 ,W2 ,W3,W4,S 
COMMON /KONS /TCRIT,PCRIT, ROCRIT, TOL, D41 R 
COHMON/ACOEF/A( 17) 

COMMON/ LCOEF/AL (17) 

COMMON/VCOEF/AV ( 17 ) 

DOUBLE PRECISION B WD , B WX , BWC VGL , BWCOND 
DOUBLE PRECISION P , T , D, ROK , DELTA , XVAL 
DOUBLE PRECISION P P , BT , C T , B P , CP , H , D1 , D2 , D3 , 04 , E X PF 
DOUBLE PRECISION CS ,CH,F , E,U,C , B , AA 

DOUBLE PRECISION AF , R F , C F , DF , E F , F F , B 1 , C 1 , C2 , CK , E 1 , EK 

DOUBLE PRECISION WKS , W 1 , W2 , W3 , W4 ♦ S 

DOUBLE PRECISION TCRIT ,PCRIT , ROCRIT ,T0L,D41R 

DOUBLE PRECISION A 

DOUBLE PRECISION AL 

DOUBLE PRECISION AV 

DATA ( AV ( I ) , 1=1,17) /82. 08 199823, 20. 62278898,-129279. 2029, 
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1 -7237230.1370, 115924274.50, -10108798750., ; 

2 317.62939700, 2581305.9670,' 241066.9065,-10703806250., 

3 1016369054000., -19384310020000., 3857308627000., 

4 -675746323600000., 14621146530000000., 52 549922 590 ., 1 800 . 100800/ \ 

5 , ( AL { I ) , 1=1,17)/ 82.08199823, 63.74020840, -353918.0407, 

6 -4810952.457, 91278833.49, -8816106422., 

7 -1283.735749, 8076213.444, 1425160.973, 6410245277., 

8108516291300., -2930340262000., -5235483345000., < 

9 -255111438000000., 47327993 10000000 . , 352232 77 740 . 0 , 1 800 . 1 00800/ 

DATA <TCRIT,PCRIT,R0CRIT,T0L )/ 32.98400, 12.7700, .0152672 DO, 

1 .00001 00/ 

END 

$ I BFTC LIQUID DECK 

SUBROUTINE LIQUID 
COMMON /L CO EF/ AL { 17) 

COMMON/ ACOEF/ A (17) 

DOUBLE PRECISION AL ,A 
DO 77 1=1,17 

77 A ( I ) = AL ( I ) 

RETURN 

END 

$ I R F T C VAPOR DECK 

SUBROUTINE VAPOR 
COMMON/ VCOh F/ A V ( 17) 

COMMON/ ACOFF/A ( 17 ) 

DOUBLE PRECISION AV,'A 
DO 77 1=1,17 
77 A ( I ) = AVI I ) 

RETURN 

END 

I BFTC PK DECK 

PRESSURE FUNCTION OF DENSITY GRAM-MOLE/CC AMD COMMON T TERMS 
DOUBLE PRECISION FUNCTION PKIRK2) 

COMMON/ACOEF/AI 17) 

COMMON/ FTS/ A F , BF ,CF , DF , E F , F F , B 1 ,C1,C2,CK,E1,EK 
DOUBLF PRECISION RK2,DEXP 
DOUBLE PRECISION A 

DOUBLF PRECISION AF , BF , C F , DF , E F , F F , B 1 , C 1 , C2 , CK , E 1 , EK 
EK= DEXP I — A C 1 7 ) * RK2*RK2 ) 

CF= C2 - EK + Cl 
EF= El * EK 

PK=U ( ( (FF*RK2+EF)*RK2 +DF)*RK2+CF)*RK2 + BF ) *RK2+ AF ) *RK2 

RETURN 

END 

I BFTC P PK DECK 

PARTIAL DERIVATIVE FOR PRESSURE AMD DENSITY 
DOUBLE PRECISION FUNCTION PPKIRK3) 

COMMQN/ACOEF/A ( 17) 

COMMON/ FTS/ A F , BF , C F , DF , EF , F F , B 1 ,C 1 , C2 ,CK , E 1 , EK 
DOUBLE PRECISION RK3 
DOUBLE PRECISION A 

DOUBLF PRECISION AF , BF , C F , DF , EF , F F , B 1 , C 1 , C2 , CK , E 1 , EK 

PPK= I C ( ( ( — 2 A(17)*FF*RK3+6.*FF)*RK3+(5.*EF -2 .*A ( 17 )*C2*EK ) ) 

1 -RK3+4. -DF ) *RK3 + 3.*CF) *RK3 + 2. -BF) * RK3 + AF 

RETURN 
END 

SIRFTC COEF DECK 

C FUNCTIONS OF TEMPERATURE FOR STATE CALCULATIONS 

SUBROUTINE COEF ( T 1 ) 

COMMON/ FTS/ A F,BF,CF,I)F,FF,FF,B1,C1,C2,CK, El, EK 
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COMMON/ ACOEF/ A (17) 

DOUBLE PRECISION AF , BF , CF * OF , EF , F F , B 1 , C 1 , C2 , CK , E 1 , 6K 
DOUBLE PRECISION A 
DOUBLE PRECISION U 
AF = A(l) * T1 

BF= AF*A ( 2 ) + A ( 3 ) + <(A(6)/T1/T1 + A(5))/T1 + A(4))/T1 
DF= A ( 9 ) * T 1 
FF= A ( 1 6 ) 

C2 = ( ( A ( 1 2 ) / T 1 + A { 1 1 ) ) / T 1 + A ( 10 ) ) / Tl/Tl 

E 1 = ( ( A ( 1 5 ) / T 1 + A { 14 ) ) / T 1 + A ( 1 3 ) )/Tl/Tl 
Cl = A ( 7 ) * AF + A(8) 

RETURN 

END 

IBFTC 8WPSAT DECK 

SATURATION PRESSURE ATM FUNCTION OF TEMPERATURE DEGREES KELVIN 
DOUBLE PRECISION FUNCTION BWPSAT ( T5 ) 

COMMON/KONS/TCRI T , PCR I T , ROC R I T , TOL ,D41R 

COMMON/ WORK S / WKS , W 1 , W2 , W3 , W4 , S 

DOUBLE PRECISION T5, ROKSL 

DOUBLF PRECISION WKS , W 1 , W2 , W3 , W4 , S 

DOUBLE PRECISION 1 C R I T , PCR I T , ROC R I T , TOL , 041 R 

BWPSAT = ROKSL (T5) 

BWPSAT = PCR I T - (((((18502401000000. * S - 1 2 34938 500000 . ) *S + 
1 34607546000.) *S - 629755380.) *S + 7312395. ) *S** ( 3)) 

RETURN 

END 

IBFTC BWROSV DECK 

SATURATED VAPOR DENSITY 
DOUBLE PRECISION FUNCTION BWR0SV(T8) 

BWROSV = BWDEN$(T8,BWPSAT(T8)-. 00001 ) 

RETURN 

END 

IBFTC ROKSL DECK 

SATURATED LIQUID DENSITY INTERNAL UNITS GR-MOLE/CC 
DOUBLE PRECISION FUNCTION ROKSL (T4) 

COMMON/ WORKS/ WKS, W1,W2,W3,W4,S 
COMMON/KONS/TCRI T, PCR IT, ROC R I T, TOL, D41R 
DOUBLF PRECISION WKS , W 1 , W2 , W3 , W4, S 
DOUBLE PRECISION TCR I T , PCR I T , ROCR I T , TOL , D4 1R 
IF ( T4.LT.TCR I T ) GO TO 7 
S = 0 . 

GO TO 8 

7 WKS- (TCRIT-T4)**(l./3.) 

S= (((-.000020693181 * WKS -.00018306903) *WKS 
1 +.0014973511) * WKS + .0062675345) #WKS 

8 ROKSL= ROCRIT + S 
RETURN 

END 

IBFTC BWTSAT DECK 

SATURATION TEMPERATURE DEGREES KELVIN 
DOUBLE PRECISION FUNCTION BWTSAT (P21) 

B W T S AT = 20 • 

W4= AL OG ( P2 1 ) 

DO 11 I T6R$= 1,10 
W3= BWTSAT + 1.0044 

DELT = ( W4 - (4.60659779 - 115.35279 / W3+ .0402605852* 

1 BWTSAT ))/( 115.35279 / W3**2 + .0402605852) 

BWTSAT = BWTSAT + DELT 
IF (DABS (DELT) - .0005) 3,3,11 
11 CONTINUE 
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3 IF ( BWTSAT *GE. 32*984) BWTSAT= 32*983 
RETURN 
END 

IBFTC DPDT DECK 

DPDT DERIVATIVE FOR THERMODYNAMIC GROUP REQUIRES BASIC MODULES 
DOUBLE PRECISION FUNCTION DPDT (T13) 

COMMON/ WORKS/WKS »W1»W2,W3,W4,S 
COMMON/ PTD/P,T,D,ROK,DELT A t XVAL 
COMMON/ FTS/AF,BF,CF, OF, EF,FF,B1 ,C1,C2,CK,E1, EK 
COMMON/ AC OEF/ A ( 17) 

DOUBLE PRECISION WKS , W 1 t W2 t W3 , W4 , S 

DOUBLF PRECISION AF , BF ,CF , DF , EF , F F , B 1 , C 1 , C2 , CK , E 1 , EK 
DOUBLE PRECISION A 

DOUBLF PRECISION P , T , D , ROK , DE L TA , X VAL 

Wl=EK*ROK 

W2=ROK**2 

DPDT = W2* (((((-4.* ( W 1 * ( A < 1 5 ) #W2 + A(I2)) + A ( 6 > ) ) 

1 / T I 3 ~ 3 * # W 1 * ( A ( 1 4 ) * W2 + A(I1)))/T13 - 2** 

2 (W1 # ( A ( 13 ) *W2 + A ( 10 ) ) + A ( 5 ) ) ) / T 1 3 

3 —A (4) )/T13/T13 + A(1)/R0K + All) * A(2) + 

4 ROK# ( A ( 1 ) *A<7) + A ( 9 ) # ROK)) 

RETURN 

END 

IBFTC BWTEMP DECK 

FUNCTION SUBROUTINE FUR TEMPERATURE IN C.G.S. UNITS 
AS FUNCTION OF DENSITY AND PRESSURE. REQUIRES BASIC GROUP 
USED FOR DENSITY CALCULATION 
DOUBLF PRECISION FUNCTION BWTEMP ( DI N, P20 ) 

COMMON/ PTD/P,T ,D,ROK, DELTA, XV AL 

COMMON/ FTS/AF,BF,CF,DF,EF,FF,B1,C1,C2,CK, El , EK 
COMMON/ACOEF/A( 17) 

COMMON/ KONS/TCR IT,PCRIT, ROCR IT,T0L,D41R 

COMMON/ BWDXCK/BWD,BWX,BWCVGL,BWCOND 

DOUBLE PRECISION B WD , B WX , BWC VGL * BWCOND 

DOUBLE PRECISION P , T , D , ROK , DE LTA , XV AL 

DOUBLF PRECISION AF , B F , C F , !)F , EF , FF , B 1 , C 1 , C2 t CK , E 1 , EK 

DOUBLE PRECISION A 

DOUBLE PRECISION TCR I T , PCR I T , ROCR I T , TOL , D4 1 R 
D20= DIN 
T — TCR I T 

ROK = D20/2. 01572 
WKS= ROK 

I F ( P20-PCR I T ) 10,13,13 
13 CALL VAPOR 
XVAL- 2. 

CALL COEF(T) 

IF (PK(ROK) - P20 ) 55,11,11 

10 T = BWTS AT ( P20 ) 

CALL VAPOR 

I F { ROK-ROCR I T ) 22,23,23 

23 ROSL= ROKSL (T) * 2.01572 
IF ( D20 - ROSL) 15,15,11 

11 CALL LIQUID 
X VA L = -1* 

GO TO 55 

22 ROS V= BWROSV (T) 

IF ( ROSV - D20) 16,16,55 
55 DO 67 1= 1,10 
ROK= WKS 
CALL COEFIT) 
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EK=DEXP(-A( 17)*RQK**2) 

DELTA= -(P20-PMR0K) ) 

DELTA= DELT A/ DPDT { T ) 

T = T— DEL T A 

I F ( DABS ( DELTA ) — .0003) 77,77,67 
67 CONTINUE 
GO TO 77 

15 ROSV = BWROSV (T) 

GO TO 100 

16 ROSL = ROKSL(T) * 2.01572 

100 XVAL = (ROSL - D20)/ (ROSL - ROSV) * R0SV/D20 
77 R WTEMP=T 
B WD= D20 
B WX= X VAL 
RETURN 
END 

IBFTC BWPRES DECK 

FUNCTION SUBROUTINE FOR PRESSURE IN C.G.S. UNITS 
PRESSURE ATM AS FUNCTION OF TEMPERATURE AND DENSITY 
REQUIRES BASIC MODULES USED FUR DENSITY CALCULATION 
DOUBLE PRECISION FUNCTION BWPR E S ( T 1 7 , DI N ) 

COMMON/ACOEF/A 

COMMON/ PTD/P,T,D,ROK, DELT A, XVAL 
COMMON/KONS/TCR IT , PCR IT ,RQCR IT , TOL , D 4 1 R 
COMMON/ BWDXCK/BWD,BWX f RWCVGL, BWCOND 
DOUBLE PRECISION BWD, B WX , BWC VGL , BWCOND 
DOUBLE PRECISION A 

DOUBLE PRECISION P , T , D , ROK , DEL TA , X V AL 
DOUBLE PRECISION TCR I T , PC R I T , ROC R I T , TOL , D4 1 R 
017= DIN 
CALL VAPOR 
X VAL= 2. 

ROK= Dl 7/2.0157 
I F ( TCR I T-T 1 7 ) 12,13,13 

13 B W P R E S = 8WPSAT ( T 1 7 ) 

I F ( ROK-ROC R I T ) 17,14,14 

14 ROSL = ROKSL (T17) * 2.01572 
IF (ROSL - D 1 7 ) 11,72,7 2 

11 CALL LIQUID 
XVAL= -1. 

GO TO 12 

17 R0SV=BWR0SV(T17) 

R0K=D17/2.0157 

I F ( ROSV— D17 )71,12, 12 

71 ROSL= ROKSL ( T 1 7 ) * 2.01572 
GO TO 100 

72 ROSV= BWROSV (T17) 

100 XVAL = (ROSL - D17) / (ROSL - ROSV)* ROSV / D17 
GO TO 77 

12 CALL COEF ( T 1 7 ) 

BWPR ES= PK(ROK) 

77 BWX=XVAL 
BWD=D1 7 
RETURN 
END 

IBFTC PARALO DECK 

PARALO SUBROUTINE FOR IDEAL COEFFICIENTS, SPECIFIC HEAT OF 

PARAHYDROGEN FOR TEMPERATURES LESS THAN 100 DEGREES KELVIN 
SUBROUTINE PARALO 
COMMON/C VS/CS,CH,F,E,U,C,B, A A 



DOUBLE PRECISION CS ,CH, F , E ,U T C , B t AA 

AA = 4*977816011000 DO 

B = 3384077523 D-2 

C = 3.521443735 D-4 

U “ -14.35633178D-6 

E = 23.03247505 D-8 

F = —10.38316229 D-10 

CH = .43 DO 

CS = -.1984 DO 

RETURN 

END 

$ I BFTC PAR AH I DECK 

C PAR AH I SUBROUTINE FOR IDEAL SPECIFIC HEAT PAR AHVDROGEN 

C COEFFICIENTS FOR TEMPERATURES ABOVE 100 DEGREES KELVIN 

C BUT LESS THAN 400 DEGREES KELVIN 

SUBROUTINE PA RAH I 
COMMON/CVS/CStCH, F,E»U*Cf B,AA 
DOUBLF PRECISION CS , CH, F , E , U , C , B v AA 
AA = -8.543812498 DO 
B = 29.35073908 D-2 
C = -19.64505888 D-4 
U = 6.145180979 D-6 
E = -.9172545570 D-8 
F = .05273284618 D-10 
CH = 432.3534 DO 
CS = 41.002 DO 
RETURN 
END 

$ I BFTC NORMLO DECK 

C NOR ML 0 SUBROUTINES FOR IDEAL SPECIFIC HEAT NORMAL HYDROGEN 

C TEMPERATURES LESS THAN 100 DEGREES KELVIN 

SUBROUTINE NORMLO 
COMMON!/CVS/CS,CH t F,EtU,C,ft, A A 
DOUBLE PRECISION CS ,CH, F t E, U,C ,B, AA 
AA = 4.960343508 DO 
B = .02245403779 D-2 
C = .4825909452 D-4 
U = -2.895260963 D-6 
E = 5.0980396 D-8 
F = -2.275937407 D-10 
CH = 254.03 DO 
CS = 4.151 DO 
RETURN 
END 

S I BFTC NORMHI DECK 

C NORMHI SUBROUTINE FOR IDEAL COEFFICIENTS SPECIFIC HEAT OF 

C NORMAL HYDROGEN TEMPERATURES ABOVE 100 DEGREES KELVIN UP TO 

C 400 DEGREES KELVIN 

SUBROUTINE NORMHI 

COMMON/C VS/CS t CH t F t E,U,C,B,AA 

DOUBLE PRECISION CS ,CH , F , E t U T C , B , A A 

AA = 3.809873570 DO 

B = 1.2237926310 D-2 

C = .8816230861 D-4 

U = -.6662777844 D-6 

E = .1533405461 D-8 

F = -.0121587109 D-10 

CH = 301.14435 DO 

CS = 8.11846885 DO 

RETURN 
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END 

$ I BFTC COEFID DECK 

C COEFID SUBROUTINE FOR SETTING UP IDEAL COEFFICIENTS 

C NORMAL OR PARA MODE 

SUBROUTINE COEFID ( T 1 ) 

COMMON/ FIXED/MODE, ITERS 
I F { MODE—1 ) 1 * 1 1 2 

1 IF(Tl-100.) 1 1 « 12 « 12 

11 CALL PARALO 
RETURN 

12 CALL PAR AH I 
RETURN 

2 1 F ( T 1—100 • ) 21*22*22 

21 CALL NOR ML 0 
RETURN 

22 CALL NORMHI 
RETURN 

END 

S I BFTC BWCPIG DECK 

C FUNCTION SUBROUTINE FOR SPECIFIC HEAT IDEAL C-AS NORMAL OR 

C PARAHYDROGEN ALL TEMPERATURES UP TO 400 DEGREES KELVIN 

DOUBLE PRECISION FUNCTION BWCPIG(TIO) 

COMMON/CVS/CS*CH*F*E*U*C*B*AA 
DOUBLE PRECISION C S * C H , F , E , U , C , B * A A 
CALL COEFID(TIO) 

BWCPIG =l./2.01572*( ( ( ( ( F * T 1 0+ E ) *T 1 0+U ) * T 1 0+C ) * T 1 0+B ) 

1 *T10+AA) 

RETURN 

END 

$ I BFTC BWHIDG DECK 

C FUNCTION SUBROUTINE FOR ENTHALPY IDEAL GAS T LESS THAN 400 DEG K 

DOUBLE PRECISION FUNCTION RWHIDG(TIO) 

COMMON/CVS/CS*CH*F,E *U*C ,B*AA 
DOUBLE PRECISION C S * C H , F , E , U f C * B , A A 
CALL COEFID (T10) 

BWHI0G=l./2.01572*<CH+( ( ( ( ( F/6. *T 10+E/ 5. ) *T1 0 
1 +U/4. ) * T 1 0+C / 3 . ) * T 1 0+ B / 2 . ) *T 1 0 + A A ) *T 1 0 ) 

RETURN 

END 

$ I B FTC BWSIDG DECK 

DOUBLE PRECISION FUNCTION BWSIDG(TIO) 

C FUNCTION SUBROUTINE FOR ENTROPY IDEAL GAS T LESS THAN 400 DFG K 

COMMON/CVS/CS*CH,F ,E *U,C * B* AA 
DOUBLE PRECISION C S * C H , F * E * U * C * B , A A 
CALL COEFID (T10) 

BWSIDG ssl./2.01572*(CS+AA*DL0G(T10)+( <<<F/5.*T10 

1 +E/4.)*T10+U/3.>*T10+C/2. ) #T 10+B ) *T 10 ) 

RETURN 

END 

I8FTC BWSSTL DECK 

FUNCTION SUBROUTINE FOR ENTROPY SATURATED LIQUID PARAHYDROGEN 
DOUBLE PRECISION FUNCTION BWSSTL(Tll) 

COMMON/KONS/TCRIT, PCRITtROCRIT, T0L,D41R 
DOUBLE PRECISION TCR I T * PCR I T * ROCR I T , TOL * D41 R 
BWSSTL = 1./2. 01572*(5. 77919623-2. *.43272654* 

1(TCRIT-T11)**(1./2.I - 1. 3333618 *DLOG( Til )+ 

2 ( ( .0001 7227638/ 3.*TU-.0 12 86952 7/2. )*T 11 

3 +.36251851 )*TU ) 

RETURN 

END 
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$ I BFTC BWHSTL DECK 

C FUNCTION SUBROUTINE FOR ENTHALPY SATURATED LIQUID PARAHYDROGEN 

DOUBLE PRECISION FUNCTION BWHSTL(Tll) 
COMMON/KONS/TCRIT,PCRIT,ROCRIT,TOL,D41R 
COMMON/ WORKS/WKS,Wl ,W2,W3,W4,S 
DOUBLF PRECISION WKS , W 1 , W2 , W3 , U4 , S 
DOUBLE' PRECISION TCR I T , PC R I T , ROCR I T t TOL , D4 1 R 
W = ROKSL(TH) 

W = DLOG(W)*3. 136424038+ 

1 S*1000.*( ( ( ( { 2 1 586 1 3. 45*S- 190 139. 7999 >*S 

2 +7954. 571 19 >*S-245. 892 7562 )*S+6. 72 800008 >*S 

3 —.2054354458 ) 

W = -120. 142080*(W+13. 073164) 

BWHSTL = W+1./2 .01572* ( -53 . 4 3794-2 •*. 432 72 654 

1 *(2.*TCRIT+Tll)/3.* (TCRIT-T11 )**( 1./2.) 

2 +( ( ( . 000 1 722 7638/4. *T11-. 01 286952 7/ 3. ) *T I L 

3 +.36251851/2. )*T1 1-1. 3333618 ) - T 1 1 ) 

RETURN 

END 

I BFTC BWSSV DECK 

FUNCTION SUBROUTINE FOP. ENTROPY SATURATED VAPOR HYDROGEN 

EITHER PARA OR NORMAL AS A FUNCTION OF SATURATION PRESSURE 
DOUBLE PRECISION FUNCTION BWSSV(P23) 

COMMON/ PTD/P,T ? 0, ROK, DELTA, XV AL 
DOUBLE PRECISION P , T , D,ROK , DELTA , XVAL 
T = BWTSAT { P23 ) 

RUK= 0 

BWSSV= SKS V ( P23 ) 

D= R0K*2. 01572 

RETURN 

END 

I BFTC SKSL DECK 

FUNCTION SATURATED LIQUID ENTROPY AS FUNCTION OF PRESSURE 
DOUBLF PRECISION FUNCTION SKSL ( P24 ) 

COMMON/ PTD/P,T,D,ROK, DELTA, XVAL 
DOUBLF PRECISION P , T t D, ROK , DELTA, XV AL 
CALL LIQUID 
SKSL= BWSSTL(T) 

RETURN 

END 

I BFTC BWSSL DECK 

ENTROPY SATURATED VAPOR INTERNAL UNITS REQUIRES BASIC AND 
THERMODYNAMIC MODULES 
DOUBLE PRECISION FUNCTION BWSSL(P25) 

COMMON/ PTD/P,T,D,ROK, DELTA, XVAL 
DOUBLF PRECISION P, T ,D, ROK, DELTA, XVAL 
T = BWTSAT ( P2 5 ) 

BWSSL= SKSL ( P2 5 ) 

D= R0KSL(T>*2. 01572 

RETURN 

END 

I BFTC HKSV DECK 

ENTHALPY SATURATED VAPOR PARA OR NORMAL HYDROGEN REQUIRES 
BASIC AND THERMO GROUP OF MODULES 
DOUBLE PRECISION FUNCTION HKSV ( P22 ) 

COMMON/ PTD/P, T,D, ROK, DELTA, XVAL 
COMMON/KONS/TCRI T , PCR I T , ROC R I T , TOL ,041R 
C 0 M M 0 N / W 0 R K S / W K S , W 1 , W 2 , W 3 , W 4 , S 
DOUBLE PRECISION P ,T ,D, ROK , DELTA , XV AL 
DOUBLE PRECISION TCR I T , PCR I T , ROCR I T , TOL , D4 1 R 
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DOUBLE PRECISION WK S T W 1 , W2 , W 3 , W4 r S 
DSV= BWROSV(T) 

12 HKSV= BWHIDG(T) + HDEV ( T , ROK , P2 2 ) 

RETURN 

END 

$ I B FTC HKSL DECK 

DOUBLE PRECISION FUNCTION HKSL(P24) 

COMMON/ PTD/P,T,D, ROK, DELTA, XVAL 
DOUBLE PRECISION P , T T l) ,ROK , DELTA , XVAL 
CALL LIQUID 
HKSL= BWHSTL ( T ) 

RETURN 

END 

SIBFTC SKS V DECK 

DOUBLE PRECISION FUNCTION SKSVIP22) 

COMMON/ PTD/P,T,D, ROK, DELTA, XVAL 

COMMON /KONS/TCR IT, PCRI T t ROCRI T,TOL ,04 1R 

COMMON/ WORKS/ WKS, W1 ,W2 , W3,U4, S 

DOUBLE PRECISION P , T , D , ROK , DEL T A , XVAL 

DOUBLF PRECISION TC R I T , PCR I T , ROCR I T , TOL , D4 1 R 

DOUBLE PRECISION WK S , W 1 , W 2 , W3 , W4 , S 

DSV= BWROSV(T) 

12 SKSV= BWSIDGIT )+SDEV( T »ROK,P22 ) 

S I DG= B W S I DG { T ) 

RETURN 

END 

IBFTC HDEV DECK 

ENTHALPY DEVIATION TERM REQUIRES THERMO AND BASIC FUNCTIONS 
DOUBLE PRECISION FUNCTION HDEV ( TO, RD, PD ) 

COMMON/FT S/ A F,BF,CF, OF, EF,FF,B1,C1,C2,CK, El ,EK 
COMMON/ ACOEF/A ( 17 ) 

DOUBLF PRECISION AF , B F , C F , DF , E F , F F , B 1 , C 1 , C 2 , C K , E 1 , EK 
DOUBLE PRECISION A 
C K= ( 1 • - EK ) / 2 • / A { 17 ) 

C 2= { CK-EK*RD*RD/2 • ) /A { 17 ) 

HOE V= 1 */2. 01572/41. 29283* ( ( < A ( 1 6 ) / 5 . *RD**3+ A ( 8 ) / 2 • > 

1 *RD+( (5.*A(6) /TD/TD+3.*A(5) )/TD+2.*A(4) )/TD+A( 3) )*RD 

2 +1./TD/TD*{ { ( 5. *A { 12 ) /TD+4.*A( 1 1 ) )/TD+3.*A{ 10 ) ) *CK 

3 +( (5.*A{ 15)/TD+4.*A( 14) )/TD+3.*A( 13) }*C2 ) 

4 +PD/RD-AF) 

RETURN 

END 

IBFTC BWHSV DECK 

DOUBLF PRECISION FUNCTION BWHSV(P23) 

COMMON/ PTD/P,T,D, ROK, DELTA, XVAL 
DOUBLF PRECISION P , T , D , ROK, DE L T A, XV AL 
T = BWTSAT(P23) 

R OK = 0 

BWHSV= HKS V ( P23 ) 

D= R0K*2. 01572 

RETURN 

END 

IBFTC BWHSL DECK 

ENTHALPY SATURATED LIQUID INTERNAL UNITS REQUIRES BASIC AND 
THERMO GROUP OF FUNCTIONS 
DOUBLE PRECISION FUNCTION F.WHSKP25) 

COMMON/ PTD/P,T,D, ROK, DELTA, XVAL 
DOUBLE PRECISION P , T , D , ROK , DE L TA , X V A L 
T = BWTS AT ( P2 5 ) 

BWHS L= HKSLIP25) 


D= R0KSL(T)*2. 01572 

RETURN 

END 

SIBFTC S DEV DECK 

C DEVIATION TERM FOR ENTROPY FUNCTION REQUIRES BASIC AND THERMO 

DOUBLE PRECISION FUNCTION SDEV ( TD, R D, PD ) 

COMMON/ FTS/ AF, BF, CF,DF, EF, FF,B1 ,C1 *C2*CK, El, EK 
COMMON/ AC OE F/ A (17) 

DOUBLE PRECISION A 

DOUBLF PRECISION AF , BF , C F , DF , E F , F F , B 1 , C 1 , C2 * CK , E 1 , EK 
CK= <l.-EK)/2./A<17> 

C2 = (CK-EK*RD*RD /2.)/A(17) 

SDEV = 1./2. 01572/41. 292 83* ( ( ( -A ( 9 ) / 3 .*RD— ( A ( 7 ) *A ( 1 ) )/2. ) 
l*RD-< <-4.*A(6)/TD/TD-2.*A(5) )/TD-A(4) ) / TO/ TD 
2-A(l )*A(2) ) *RD-A ( 1 ) *DLOG ( RD*A < 1 )*TD)+ 

3 ( ( (4.*A< 12)/TD+3.*A( 11 ) ) / TD+2 . *A ( 1 0 > )*CK+ 

4 { (4.*A( 15)/TD+3.*A( 14)) /TD+2.*A( 13) ) *C2 ) / TD/ TD/ TD ) 

RETURN 

END 

SIBFTC BWENHY DECK 

C ENTHALPY OF PARA OR NORMAL HYDROGEN AS FUNCTION F TEMPERATURE 

C AND PRESSURE IN C.G.S. UNITS REQUIRES BASIC AND THERMO GROUPS 

DOUBLE PRECISION FUNCTION BWENH Y ( T 1 2 , P 12 ) 

COMMON/ PTD/P, T,D,RDK, DELTA, XV AL 
COMMON/ WORKS/WKS, W1 ,W2,W3,W4,S 
CQMMON/ACOEF/A (17) 

COMMON/ FTS/ AF,BF,CF,DF,EF,FF,B1,C1,C2,CK, E1,EK 
DOUBLF PRECISION AF , BF , CF , DF , E F , F F , B 1 , C 1 * C2 , CK , E 1 , EK 
DOUBLE PRECISION A 

DOUBLF PRECISION P , T , D , ROK , DEL TA , X VAL 
DOUBLE PRECISION WK S , W 1 T W2 , W3 , W4 , S 
D= B WDENS ( T1 2 * P 1 2 ) 

BWENHY = HDEV(T12,ROK,P12 ) 

IF(XVAL-2.) 1*2,2 
2 BWENHY=BWENHY+ BWH I DG ( T 12 ) 

RETURN 

1 W4= ROKSL ( T 1 2 ) 

EK= DEXP (-A(17)*W4*W4) 

W3= HDEVCT12, W4 , BWPSAT(T12) ) 

PSA=BWPSAT ( T12 ) 

HL=BWHSTL ( T 12 ) 

BWENHY= BWENHY+BWHSTL(T12)-W3 

RETURN 

END 

SIBFTC BWERPY DECK 

C ENTROPY OF PARA OR NORMAL HYDROGEN AS FUNCTION F TEMPERATURE 

C AND PRESSURE IN C.G.S. UNITS REQUIRES BASIC AND THERMO GROUPS 

C REQUIRES THERMODYNAMIC AND BASIC MODULES 

DOUBLE PRECISION FUNCTION BWE RP Y ( T 1 2 , P 1 2 ) 

COMMON/ PTD/P* T*0, ROK, DELT A, X VAL 
COMMON/ WORKS /WKS,W1 ,W2*W3,W4,S 
COMMON/ ACOEF/ A ( 17) 

COMMON/ FTS/ AF,BF,CF*DF,EF,FF,B1 ,C1 ,C2,CK, E1 T EK 
DOUBLE PRECISION AF , BF * C F , DF , E F * FF , B 1 * C 1 , C2 * CK , E 1 , EK 
DOUBLE PRECISION A 

DOUBLE PRECISION P , T , D , ROK , DEL T A , X V AL 
DOUBLE PRECISION WKS * W 1 , W2 * W3 , W4 , S 
D= BWDENS ( T 1 2 * P 1 2 ) 

B W E R P Y = SDEV(T12*R0K,P12) 

I F ( XVAL-2 . ) 1*2,2 
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2 8WERPY= BWERPY+BWSIDG(T12) 

RETURN 

1 W4 = ROKSL ( T 12 ) 

EK = DEXP <-A(17)*W4*W4) 

W3= SDEV(T12,W4,BWPSAT(T12) ) 

BWERPY= BWERPY+BWSSTL(T12)-W3 

RETURN 

END 

IBFTC BWTOFH DECK 

SUBROUTINE FOR TEMPERATURE AS FUNCTION OF ENTHALPY AND PRESSURE 
IN C .G.S » UNITS FOR NORMAL OR PARAHYDROGEN REQUIRES BASIC 
AND THERMODYNAMIC MOOULFS 
DOUBLE PRECISION FUNCTION BWTO FH { H2 3 , P2 3 ) 

COMMON/ B W DXCK/ B WD, 8WX, B WC V GL ,BWCOND 
COMMON/ PTD/ P,T,D,ROK, DELTA r XVAL 
COMMON/KONS/TCRIT,PCRIT,ROCR IT,T0L,D41R 
DOUBLE PRECISION BWD , BW X , 8WCVGL , BWCOND 
DOUBLE PRECISION P , T T D, ROK t OELTA , XVAL 
DOUBLE PRECISION TC R I T , PCR I T , ROCR I T , TOL , D41 R 
IF(P23-PCRIT>1,1,2 

1 T =BWTSAT ( P23 ) 

R0K = 0 

HSV=HKSV( P23) 

HSL=HKSL ( P23 ) 

IF ( HSL - H23) 4,3,3 

3 H=HSL 
GOTO 11 

4 IF ( H2 3 - HSV) 8,8,9 

8 BWTOFH = T 

BWX= (H23 - HSL) / (HSV - HSL) 

D= ROKSL ( T ) 

BWD= 2.01572* (ROK* D) / ( BWX* D + (1. - 6WX ) *ROK ) 

RETURN 

2 T = 33.0 

H=BWENHY(33.0,P23) 

IF (H23 - H) 29,29,19 

9 H = HSV 
GOTO 19 

29 T -32 • 9 

H = ft W E N H Y ( T , P2 3 ) 

IF (H23 - H) 11,30,30 

30 BWTOFH = TCRIT 
RETURN 

19 TNU = T 

TNI = TN0+10. 

R EG= 2. 

TSP = TNO 
GO TO 49 
11 TNU = T 

TNI = TNO - 4.0 
REG= -1. 

TSP= TNO 

49 FT 0 = H - H23 

DO 77 I T E R = 1,10 

FT 1 = BWENHY (TN1,P23) - H23 

DELTA = ((TNI - TNO)* FT 1 ) / ( FT 1 - FTO ) 

TN2 = TNI - DELTA 

IF ( DABS ( FTO ) - DABS ( FT1 ) ) 71,72,72 
72 TNO = TNI 
FTO = FT 1 
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71 TNI = TN2 

DELTA=DABS(TN2-TN0) 

IF (DELTA - .001 ) 76,77,77 
77 CONTINUE 
76 BWTOFH = TN2 

IF ( ABS(REG-BWX) .LT. .001) GO TO 277 
IF (REG.LT.2.) GO TO 177 
BWD= BWROSV(TSP) 

178 BWX = REG 

BWTOFH= TSP 
277 RETURN 

177 BWD= 2.01572 s !* ROKSL(TSP) 

GO TO 178 
END 

IBFTC BWTOFS DECK 

SUBROUTINE FUNCTION FOR TEMPERATURE GIVEN ENTROPY AND PRESSURE 
FOR NORMAL /OR PARAHYDROGEN IN C.G.S. UNITS REQUIRES BASIC AND 
THERMODYNAMIC FUNCTIONS 
DOUBLE PRECISION FUNCTION BWTOFS I S23, P23 ) 

COMMON/ BWDXCK/BWD,BWX,BWCVGL,BWCOND 
COMMON/ PT D/P,T,D,ROK, DEL T A, XVAL 
COMMON/KON$/TCRIT,PCRIT,ROCRIT,TOL,D41R 
DOUBLE PRECISION BWD , BW X , BWC VGL , BWCOND 
DOUBLE PRECISION P , T , D, ROK , DEL T A , X V AL 
DOUBLE PRECISION TCR I T , PCR I T , ROC R I T , TOL , D41R 
IF( P2 3- PC R I T ) 1,1,2 

1 T = BWTSAT ( P2 3 ) 

ROK- 0. 

SSV= SKS V { P2 3 ) 

SSL= SKSL ( P2 3 ) 

I F ( SSL-S23 ) 4,3,3 

3 S= SSL 
GO TO 11 

4 IF(S23-SSV) 8,8,9 

8 BWTOFS= T 

BWX = (S23 - SSL) / (SSV-SSL) 

D= ROKSL(T) 

B WD= 2.01572* (ROK* D) / (BWX* D + (1. - BWX ) * R 0 K ) 

RETURN 

2 T = 33 • 0 

S=BWERPY(33.0,P23) 

I F ( S23-S ) 29,29,19 

9 S= SSV 
GO TO 19 

29 S=BWERPY(32.5,P23) 

I F ( S23-S ) 39,30,30 

30 BWTOFS= TCR I T 
RETURN 

39 T = 32 . 5 
GO TO 11 
19 TNQ= T 

TN 1 = TN0+10. 

REG= 2 . 

TSP= TNO 
GO TO 49 
11 TNO= T 

TNI = TNO— 4 .0 
REG= -1. 

TSP= TNO 
49 FTO= S-S23 
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DO 77 I TER= 1,10 

FT 1 = BWERPY(TN1 ,P23)-S23 

DEL T A = { (TN1-TN0)*FT1)/(FT1-FT0) 

TM2= TNI -DEL T A 

IF(DABS(FTG) -DABS ( FT1 ) ) 71,72,72 
72 T N0= TNI 
FTO= FT 1 
71 TN 1 = TN2 

DELTA=DARS { TN2-TN0) 

IF (DELTA - ,001 ) 76,77,77 

77 CONTINUE 
76 BWTOFS= TN2 

IF ( ABS ( REG-BWX ) .LT. .001) GO TO 277 
IF (REG • L T • 2 • ) GO TO 177 
B WD= BWROSV ( TSP ) 

178 BWX = REG 

B WTDFS= TSP 
277 RETURN 

177 BWD=* 2 .01572# ROKSL(TSP) 

GO TO 178 
END 

$ I B FTC BWFD2 DECK 

C SECOND-DERIVATIVE FUNCTION FUR SPECIFIC HEAT CALCULATION 

C REQUIRES BASIC AND THERMODYN AN I C MODULES 

DOUBLE PRECISION FIJNC1 ION RWFD2 ( T24 , R02 A ) 

COMMON/ACOFF/A ( 17 ) 

COMMON/ FTS/AF,BF,CF,DF,FF,FF,B1 ,C1 »C2 ,CK , El , EK 
DOUBLE PRECISION A 

DOUBLE PRECISION AF , BF , CF , DF , EF , FF , 3 1 ,C 1 ,C2 ,CK , E 1 , EK 
EK= DEXP(-A(17)*R024**2) 

BWFD2= 1 •/T24**3/ 83.238* ( ( 10 **A ( 6 ) / T24**2+3 .* A ( 5) 

1 +A(4)*T24)*2.*R024+ 

2 ( ( 1 .-EK ) / A( 17)*( ( 1().*A( 12 >/T24+6,*A( 11 ) ) / T24 + 3 . *A ( 10) ) )+ 

3 ( ( l.-(A(17)*R024**2+l. ) *EK ) / A ( 17 ) **7) * 

A ( ( 10.# A ( 15)/T24+6.*A( 14) ) / T24+3 .*A ( 1 3 ) ) ) 

RETURN 

END 

SIBFTC BWCPGL DECK 

C SUBROUTINE FUNCTION FOR SPFCIFIC HEAT CP ANL) CV 

C ALSO CALCULATES SONIC VELOCITY CO FOR PARA OR NORMAL 

C HYDROGEN IN C.G.S. UNITS REQUIRES BASIC AND THERMO MODULES 

DOUBLE PRECISION FUNCTION BWCPGL ( T2S, P25 ) 

COMMON/ BWSONV/RWCO 

COMMON/ BWDXCK/BWD,BWX, BWCVGL , BWCDND 
COMMON/ P T D/P , T ,D,ROK, DELTA, XV AL 
DOUBLE PRECISION BWCO 

DOUBLE PRECISION BWD , BWX , BWCVGL , BWCOND 
DOUBLE PRECISION P, T , D , ROK , DELTA , XVAL 
0= BWDENS(T25,P25) 

I F ( XVAL — 2 • ) 11,22,22 

22 IF (T25.LE.46.) GO TO 10 

23 BWCVGL- BWCP I G ( T2 5 ) - .98584 -BWF D2 ( T2 5 , ROK ) 

GO TO 66 

10 IF (D.LT. .030774) GO TO 23 

11 R WC VGL = BWFD2 ( T25, .03782 1 ) 

B WC VGL= BWCVGL— BWFD2 ( T25 , ROK ) 

I F ( 42 • 5— T2 5 ) 3,12,12 

12 I F ( 24 . 5— T 2 5 ) 2,1,1 

1 B WCVGL=ft WC VGL+ ( ( - . 0 1 1 *T2 5+0 . 6 7 ) * T2 5+2 . 32 7 ) / 8 . 4337 
GO TO 66 
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2 BWCVGL= BWCVGL+( ( -.003*T2 5+0. 302 ) *T2 5+6 . 554 ) / 8 . 4337 
GO TO 66 

3 BWCVGL=BWCVGL+( ( .0002 *125+. 0492 ) *T25+1 1 .539 8 ) / 8. 43 37 
66 WKS= DABS ( PPK(ROK) ) 

BWCPGL= BWCVGL + T25 * 0PDKT25) **2 /ROK**2 / WKS / 83.238 
BWC(J= DABS { WKS * BWCPGL / bWCVGL *502673.98) 

BWCO= SORT ( BWCO ) 

RETURN 

END 

$ I BFTC BWVISC DECK 

VISCOSITY SUBROUTINE AS FUNCTION OF TEMPERATURE AND PRESSURE 
FOR NORMAL OR PARAHYDROGEN ALSU CALCULATES CONDUCTIVITY K 
REQUIRES BASIC AND THERMO GROUPS AND ALSO HIGH TEMPERATURE 
DOUBLE PRECISION FUNCTION B W V I SC ( T 1 5 , P I 5 ) 

COMMON/PTD/P, TtDtROK » DELTA, XV AL 
COMMON/ RWDXCK/B WD, BWX , BWCVGL , SWCOND 
DOUBLE PRECISION P , T ,D,ROK , DEL r A , XVAL 
DOUBLE PRECISION BWD, BWX, BWCVGL , BWCOND 
T 400= 400. 

I F ( T 1 5.GT. T400 ) GO TO 401 
CPIG=BWCPIG(T15> 

50 D= BWDENS ( T15, PI 5 ) 

GO TO 51 

401 CP I G= HTCPIG(T15) 

D= HTDENS ( T15tP15) 

51 CONTINUE 

VI =0.8411* l.D-4 * ( . 1 0 1 7* T 1 5 ** ( 3 . / 2 • ) * 

1 ( T 1 5 + 650.39) )/( T15 + 1 9 . 5 5 ) / ( T 1 5+ 1 1 7 5 . 9 ) 

OVTERM= -58.75 * (D/.07)**3 
IF ( OVTERM. LT. -67. ) OVTERM = -67. 

DTERM= 5.7694 + 65. * D**<3./2.) - . 000006*DEX P ( 1 2 7 . 2*D ) 

1 + 1./T15* (10. + 7.2 * ( { D / • 0 7 ) * * 6 

2 -<D/.07)**(3./2.) ) - 17.63 * DEXP (OVTERM)) 

BW V I SC= 1.D4 * (VI + 1.0-6 *D *DEXP (DTERM)) 

BWCOND= 1.D4 *VI*( (-T15*. 004458+ 1 . 8 34 1 } / 2 . 0 1 572 + 

1 ( .0008973*T15 + 1 . 1 308 ) *CP I G ) / ( 1 . +3. 2/T 1 5 ) + 

2 1 .04* ((( ((( ( (33615000. *D - 11243300. )*D + 1567680. )*D- 

3 116927. )*D + 49 52 *28) *D - 1 1 5 . 024 ) *D+ 1 . 22648 ) *D + 

4 .001102 )*0+ .00000184) 

RETURN 
END 

S I BFTC HTCPIG DECK __ 

HIGH TEMPERATURE SPECIFIC HEAT FUNCTION 400 - 1500 DEG K 
PARA OR NORMAL HYDROGEN 
FUNCTION HTC P I G ( T3 ) 

COMMON/ F I XED/MUDE, ITERS 
I F ( 1— MODE ) 1,2,2 
IF(700.1-T3) 2,3,3 

HTCP I G= ( .000000 107 l*T3-.00002 21 )*T 3+3. 4528 
RETURN 

I F ( T 3—700 • > 4,4,5 

HTC P I G= 1 ./2. 01 572* ( ( . 00000 1*T3-. 001005 ) *T3 + 7. 25 1 2 5 ) 

RETURN 

HTCP I G= 1 ./2. 01 572- ( ( . 0000005*T3- . 0002458 ) *T 3 + 6. 9648 1 ) 

RETURN 
END 

$ I BFTC H1HIDG DECK 

C HIGH TEMPERATURE ENTHALPY IDEAL GAS ABOVE 400 DEG K 

FUNCTION HTH I DG ( T3 ) 

COMMON/ FIXED/ MODE, ITFRS 


C 

C 


1 

3 

2 

4 

5 
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I F ( T 3—700 « > 11,11,21 

11 HTHI 0G= l./2.01572*( ( ( .000001/3.^3—.. 001005/2 • ) *T3 
1 +7.25125)*T3-112.2412> 

RETURN 

21 HTH I DG= l./2.01572*( (< .0000005/3. *T3-. 0002458/ 2 .) *T3 
1 +6.96481 )*T3-40. 5934) 

RETURN 

END 

IBFTC HTSIDG DECK 

HIGH TEMPERATURE ENTROPY IDEAL GAS ABOVE 400 DEG K 
FUNCTION HTS I DG ( T3 ) 

COMMON/ FIXED/ MODE, ITFRS 

COMMON/ BWDXCK/BWD , B WX , BWC VGL , BWCOMD 

COMMON/ PTD/P , T,D,R OK, DELTA, XVAL 

COMMON/ P VS /PP,BT ,C T , B P , C P , H , 0 1 , D2 , D3 , D4 , EX PF 

COMMON/ WORKS/ WKS,W1,W2,W3,W4,S 

COMMON /KONS/TCR IT ,PCRIT ,ROCRIT ,T0L,D41R 

DOUBLE PRECISION TCR I T , PCR I T , ROCR I T , TOL , D41 R 

IF(l-MOQE) 1,1,2 

1 HT S I DG= -4.903 + 2.76/2.01572 
GO TO 3 

2 HTS I DG= -4.903 

3 IF ( T 3 — 700 • ) 4, A, 5 

4 HT S I DG= HTSIDG+ 1 . / 2 . 0 1 572* ( 7 . 25 1 2 5*AL OG ( T3 ) 

1 +{ .000001 /2„*T3— .001005 )*T3) 

GO TO 6 

5 HTS I DG= HTSIDG+.680+ 1 . / 2 . 0 1 5 72* ( 6 . 9648 1 *ALOG ( T 3 ) + ( .0000005/ 2 . 

1 *T3-.0002458)*T3) 

6 RETURN 
END 

IBFTC HTHDEV DECK 

THERMAL DEVIATION TERM FOR ENTHALPY FUNCTIUN ABOVE 400 DEG K 
FUNCTION HTHDEVI TD,RD,PD) 

C OMMON / F I X E D/ MODE , I T E R S 

COMMON/ BWDXCK/BWD , BWX , RWCVGL , B WOUND 

COMMON/ PTD/P, T,D,R OK, DELTA, XVAL 

COMMON/ PVS/PP,BT ,CT ,BP,CP,H,D1,D2,D3, D4 , EXPF 

COMMON/ WORKS/ WKS, W1 ,W2,W3, W4,S 

COMMON /KONS/TCR I T , PCR I T , ROCR I T , TOL , 1)4 1 R 

DOUBLE PRECISION TCR I T , PCR I T , ROCR I T , TOL , D41 R 

BP = W4/TD*( { 1 .25*.22004*W2+.75*.036877)*W2-.25*.0055478 )*W3 

C P= W2/ TD**2*W2*( 2. *.04053*W2- 1.5*. 004788 )*W3**2 

F 1 = BT **2 

Dl= BT 

D2= F1/2.+CT 

D3= BT*Fl/6.+BT*CT 

D4 = Fl**2/24.+Fl*CT/2.+(CT**2)/2 . 

F 1 = BP 

F 2 = ( B P#D 1 +C P ) / 2 . 

F3= ( BP*D2+C P^Dl ) / 3 . 

F4= (BP*D3+CP-02 1/4. 

F5= <BP*D4+CP*D3)/5. 

F6= (CP*D4)/6. 

H T H D E V = 1./41 .292 8 3*1-82.082/ 2.01 572*TD*TD* 

1 ( ( ( ( ( (F6*RQK+F5>*R0*+F4)*R0K+F3)*R0K 

2 +F2)*R0K+F1 )*ROK ) +Pf>/ ROK-82 .082/2. 01 572* TD) 

RETURN 

END 

IBFTC HTENHY DECK 

ENTHALPY OF PARA UR NORMAL HYDROGEN TEMPERATURES ABOVE 400 DEG K 
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C REQUIRES HIGH TEMPERATURES MODULES 

FUNCTION HTENHY(T2,P2) 

COMMON/ FIXED/MODE, ITERS 
COMMON/ BWDXCK/BWD,BWX,BWCVGL,BWCOND 
COMMON/ PTD/P,T,D,ROK, DELTA, XVAL 
COMMON/ P VS/ PP,BT,CT,RP,CP,H,D1 ,02,03,04, EXPF 
COMMON/ WORKS/ WKS,W1 ,W2,W3,W4,S 
COMMON/KOIMS/TCRIT,PCRIT,ROCR it,tol,04ir 
DOUBLE PRECISION TCR I T , PCR I T , R0CR I T , TOL , D4 1 R 
n= HTDEMSI T2,P2) 

HTENHY= HTHDEVt T2»ROK,P2 )+HTHIDG(T2 ) 

RETURN 

END 

STBFTC HTSDF V DECK 

C ENTROPY THERMAL DEVIATION FOR HIGH TEMPERATURE FUNCTION 

FUNCTION HTSDEV(TD,RD,PO) 

COMMON/ FIXED/MODE, I TERS 

COMMON/ BWDXCK/BWD, BWX,BWCVGL,BWCOND 

COMMON/ WORKS/ WKS, W1 ,W2,W3,W4,S 

COMMON/ PTD/P,T, 0,R0K, DELTA, XV AL 

COMMON/ PVS/PP,BT , C T , BP ,C P , H , D 1 , 02 , 03 , D4 , EX PF 

C OMMON/ KONS / TCR IT, PC RIT, ROC PIT , TOL » D41R 

DOUBLE PRECISION TCR I T , PC R I T , ROCR I T , TOL , D4 1R 

8 P = W4/TD*( ( 1.25*.22004*W2+.7S*.036877 ) *W2-. 2 5*. 00f>547« ) *W3 

CP = W2 / TD**2 *W2* ( 2. *. 040 B3*W?-1 .5*. 00478 8 ) *W 3**2 

F 1 = BT**2 

01= BT 

D2 = F 1 / 2 ♦ +C T 

D3 = BT*Fl/6.+BT*CT 

04= Fl**2/24.+Fl*CT/2.+(CT**2)/2 • 

Fl= BP 

F?= ( BP*Dl+CP>/2. 

F3= ( BP*D?+C P*D1 ) / 3 . 

F4 = (BP*D3+CP*D2)/4. 

F5= ( B P*D4+CP*D3 ) / !3 • 

F6= ( C P*D4 ) / 6 • 

HTSDEV= -D41R*( ( [ ( 04/ 4 . *R0K+D3/ 3 * ) SR0K + D2 / 2 • ) *ROK 

1 +D1)*R0K 

2 +TD*( ( ( ( ( IF6#RUK+F5)*R0K+F4)*R0K+F3)*RUK 

3 + F 2 ) * R 0 K + F 1 )*ROK) 

4+ALnG<R0K*041R*r0) ) /41 .29233 

RETURN 

END 

$ I B F TC HTERPY DECK 

C ENTROPY OF PARA OR NORMAL HYDROGEN FOR CEMPERATIJRFS ABOVE 

C 400 DEG K REQUIRES HIGH TEMPFRATIJRF MUDDLES 

FUNCTION HTERPY (T2,P2) 

COMMON/ F I XED/MODF, ITERS 
C OMMON / BWDXCK/BW!) ,B W X , BN C VG L , B W C I IN D 
COMMON/ PTD/P,T,D t ROK, DELTA, XV AL 
COMMON/ P VS /PP,BT , C T , B P , C P , H , D 1 , 02 , 03 , D4 , E X P F 
COMMON/WORKS/WKS, W 1 , W2 , W3 , W4 , S 
COMMON/KOMS/TCRI T, PCR IT, ROCR II ,TllL ,D41R 
DOURLF PRECISION TCRIT ,PCRIT,ROCRIT,TOL, 04 1 R 
0= HTDFNS(T2,P2) 

H T E R P Y = HTSDEV(T2,R0K,P2)+HTSI0G( T2 ) 

RETURN 

FND 

SIBFTC HT10FS DECK 

C TEMPERATURE AS FUNCTION OF ENTHALPY AND PRFSSURE 
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C FOR PARA OR NORMAL HYDROGEN IN REGION ABOVE 400 DEG K 

FUNCTION HTT OFS ( S5 , P 5 ) 

COMMON/ P VS /PP, BT r CT ,BP,CP,H,D1 , D2 , D3 , 04 , FX PF 
COMMON/ PTD/P, T, D,ROK, DELTA, XV AL 
T = 400 • 

DO 77 1=1,10 
S=HTERPY(T,P5) 

CVGL=HTCPIG<T)~.9854- 1.98 646/82. OP 2*HTFD2 ( T , ROK ) 

D= (CVGL+T*(HTDPDT(T) /ROK ) **2/ PP/R2 . 082 1*1* 98646)/ T 
DELTA* (S-S5)/D 
T =T-DELTA 

IF( A3S( DELTA)-. 01 ) 2,2,77 
77 CONTINUE 
2 HTT OFS = T 
RETURN 
END 

MBFTC HTDPDT DECK 

C DERIVATIVE DPDT FOR HIGH TEMPERATURE REGION FUNCTIONS 

FUNCTION HTDPDT ( TD ) 

COMMON/ F IXED/MODE, ITERS 
COMMON/ BWQXCK / R WO ,BWX , BWCVGL , BWCOND 
COMMON/ PTD/P, T,D,ROK, DELTA, XV AL 
COMMON/ PVS/PP, RT , C T , B P , C P , H , D 1 , 02 , D3 , D4 , FXPF 
COMMON / WOR K S / WK S , W 1 , W? , W 3 , W4 , S 
COMMON/KONS / TCRIT ,PCRIT , ROC KIT , TOL , D4 1 R 
DOUBLE PRECISION TCR I T , PCR I T , ROCR I T , TOL , D4 1 R 
R P = 1 • / TD*W4* ( ( 1 .25*.2?004*W2+. 75*. 036877)*W2 
1 -.25*. 0055478 )*W3 

CP= 1 ./TD**2*W2*< 2. *.04053*W2-1 .5*. 004788 )*W3**2 
HTDPDT* D41R*R0K* ( ROK*TD* ( 8P+C P*ROK ) + 1 , ) *EXPF 
RETURN 
END 

SIRFTC HTDENS DECK 

C DENSITY FUNCTION GIVEN TFMPFRATURE AMD PRF.SSURE IN C,G.S. UNITS 

C FOR HIGH TEMPERATURE PARA OR NORMAL HYDROGEN 

FUNCTION HTDENS < T 1 , P 1 ) 

COMMON/ FIXED/MODE, I TFRS 

C OM MOM /BWOXCK/BWD, RWX, BWCVGL , BWCOND 

COMMON/ PTD/P, T ,D,R OK, DELTA, XVAL 

COMMON/ P VS/ PP,BT,CT, BP, C P, H, D1 ,D2,D3,D4, EXPF 

COMMON/ WORKS /WKS,W1 ,W2 ,W3,W4,S 

COMMON/ KONS/ TCR IT, PCR I T , ROCR I T , TOL , D41R 

DOUBLE PRECISION BW D , 3W X , BWCVGL , BWCOND 

DOUBLE PRECISION TCR I T , PC R I T , ROCR I T , TOL , D4 1 R 

041 R= 82.082/2.01572 

ROK= P1/T1/D41R 

W2=1./T1**( 1 . / 2 • ) 

W4= W2**< 1./2. I 
W3= 1000./ .089888 

BT= W4*( . 00 55478+ W2*(-.22004*W2~. 036877) )*W3 
CT= 1 ./T1*W2*(-.04053*W2+.0047880)*W3**2 
DO 1 1=1,10 

EXPF= EXP( (CT*R0K+BT)*R0K) 

P= D41R*T1*R0K*EXPF 

PP= D41R*T1*EXPF*(R0K*(BT+2.*CT*R0K>+1. ) 

DELT A= (P-PD/PP 
ROK= ROK “DELTA 

IF< ABS(DELTA)-. 0000005) 70,70,1 
1 CONTINUE' 

70 D=R0K 
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HTDENS= D 
BW0= D 
BWX = 2. 

RETURN 

END 

IBFTC HTFD2 DECK 

SECOND DERIVATIVE FUNCTION FOR SPECIFIC HEAT CALCULATION 
FUNCTION HTF02 (T9 * R9 } 

COMMON/FIXED/MODE, ITERS 
COMMON/ BWDXCK/BWD,BWX,BWCVGL,BWCOND 
COMMON/ PTD/P,T,D,RUK , DELTA, XV AL 
COMMON/ P VS/ PP,BT,CT,BP,CP,H,D1, 02, D3,D4,EXPF 
COMMON / W ORK S / WK S , W 1 , W 2 , W 3 , W 4 , S 
COMMON/ KONS/ TCR I T , PC R I T , ROCR I T , TOL , D4 1 R 
DOUBLE PRECISION TCR I T , PCR I T , ROCR I T , TOL , D4 1 R 
BPP= l./16./T9/T9*W4*( ( -45 . * . 22004*W2 
I -21 .*.036877 )*W2 + 5.*. 0055478 )*W 3 
CPP= 1 • /T9**3* ( — 6 • *• 04053*W2+3. 75*. 004788 ) *W2*W3*W3 
A0= 2.*BP+T9*BPP 
A 1 — BP**2+2.*CP+T9*CPP 
A 2 = 2.*T9*BP*CP 
A3= T9*CP**2 
Gl= AO 

G2= ( A0*D1+A1 )/2. 

G3= { A0*D2+A1*D1+A2 )/3. 

G4= ( A0*D3+A1*D2+A2*D1+A3 )/4. 

G5 = ( A0*D4+A1*U3+A2*D2+A3*D1 ) / 5 - 
G6= ( A1*D4+A2*D3+A3*D2 )/6. 

G7= (A2*D4+A3*U3)/7. 

G8= { A3*D4 ) / 8 • 

H T F D 2 = 82.082/2.01572*T9*< {((((( ( G8*R9+G7 ) *R9+G6 ) 

1 *R9+G5)*R9+G4)*R9+G3)*R9+G2 )*P9+G1 )*R9) 

RETURN 

END 

IBFTC HTTOFH DECK 

TEMPERATURE AS A FUNCTION OF ENTHALPY AMD PRESSURE FOR 
PARA OR NORMAL HYDROGEN ABOVE 400 DEGREES KELVIN 
FUNCTION HTTOFH(H5,P5) 

COMMON/ PTD/P ,T,0, RUK, DELTA, XV AL 

COMMON/ P VS / P P,BT,CT, BP, CP,H,D1 , D2 , D3 , D4 , E X p F 

T = 400. 

DU 77 1= 1,10 
H=HTENHY(T,P5) 

3WCVGL= HTCPIG (T)-. 9854- 1.98646/ 8?.082*HTFD2 (T,ROK) 

D= BWCVGL+T*(HTDPDT ( T ) /ROK ) **2/ P P/8 2 .0821*1 .98646 
DELTA= ( H-H5 ) / D 
T = T-DELTA 

IF(ABS( DELTA)— .01 ) 2,2,77 
77 CONTINUE 
2 HTTOFH= T 
RETURN 
END 

BFTC HTCPGL DECK 

SPECIFIC HEAT OF NORMAL OR PARA HYDROGEN AS FUNCTION 
OF TEMPERATURE AND PRESSURE IN C.G.S. UNITS ABOVE 400 DEG K 
FUNCTION HTCPGL ( T4,P4) 

COMMON/ FI XE IV MODE, ITERS 
C 0 M M ON / B W S ON V / B WC 0 

COMMON/ BWDXC K / BWD ,BWX , BWCVGL » BWCOND 
COMMON/ PT D/ P, T,D, ROK, DELTA, XV AL 
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C0MM0N/PVS/PP,BT,CT,BP,CP,H,D1,D2,D3,D4,EXPF 
COMMON/ WORKS/ WKS,W1,W2,W3,W4,S 
C0MM0N/K0NS/TCRIT,PCRIT,R0CRIT,T0L,D41R 
DOUBLE PRECISION BWCO 

OQIJHL E PRECISION BWD, BWX ,BWCVGL ,BWCOND 
DOUBLE PRECISION TCR I T , PCR I T , ROCR I T , TOL , 04 I R 
0= HTDENS ( T4 , P4 ) 

BWC VGL - HTCPIGt T4 ) -.9854-1 . 98646/ 82 . 082*HTFD2 ( T4, ROK ) 

HTCPGL = BWCVGL+T4*(HTDP0T (T4)/R0K)**2 
I /PP/82. 082*1. 9864ft 

BWCO = SORT ( PP * HTCPGL/BWCVGL * 502673.98* 2.01572) 

RETURN 

END 

SI3FTC HTTMP DECK 

FUNCTION HTTEMP(DIN r P20) 

COMMON/ BWDXCK/R WO, BWX , BWCVGL , BWCOND 
DOUBLE PRECISION BWD , BW X , BWCVGL , BWCOND 
COMMON/ PVS/PP,BT,CT,BP,CP,H t Dl , D2 , D3 , D4, EX PF 
C OMMON/ PTD/P,T,D, ROK, DELTA, XV AL 
R OK = D I N 
T=^00. 

X VAL=2 . 

DO 67 1=1,10 

DELTA = -<P20 - HT PR F S ( T , ROK ) ) 

DELTA = DELTA / HTDPDT(T) 

T = T - DELTA 

IF (ABS( DELTA) - .0003 ) 77,77,67 
67 CONTINUE 
GO TO 77 
77 HTTEMP = T 
R Wl)=D I N 
BWX= XV AL 
RETURN 
END 

SIBFTC HTPKS DECK 

FUNCTION HTPRES ( T 1 , DO ) 

COMMON/ PTD/P,T,D, ROK , DELTA, XV AL 

COMMON/ P VS/ PP,BT,CT,BP,CP,H, 01 , D2 , 03 , DA , E X PF 

D41R = 82.082/2.01572 

ROK = DO 

W2=S0RT ( 1 . / T 1 ) 

W 4 = S 0 R T { W 2 ) 

W3=1000./. 089888 

B T = W 4 * ( .0055478 +W2 *(-.22004*W2 - .036877))*W3 
C T= 1 ./T1*W2*( — .04053 *W2 + . 004788 ) *W 3**2 
E X P F = FXP( ( C T * R 0 K + BT ) *ROK ) 

UTPRES= D41R * T1 * ROK * EXPF 

RETURN 

END 

IBFTC BT DECK 

FUNCTION SUBROUTINE REQUIRES FOR BRITISH THERMAL UNITS 
PARA OR NORMAL COMPOSITION ALL PROPERTIES N = 1 - 22 PARA 
AND N = 201 - 222 NORMAL HYDROGEN. PRESSURES PSI, AMD 
TEMPERATURES DEGREES RANKIN 
FUNCTION BT(N,P1,P2) 

COMMON/FIXED/MODE, ITERS 
COMMON/ B NS ONV/ BWCO 

COMMON/ BTDXCK/BTD, BTX,BTCVGL,BTCOND,BTCO 
C OMMON/ R W DXC K / B WD , B W X , B WC VGL , BWC OND 
DOUBLE PRECISION BTD , BT X , BTC VGL , B TCOND , B TC 0 



DOUBLE PRECISION BWD , BW X , BWCVGL , BWCOND 

DOUBLE PRECISION BWCO 

MODE = N/lOO 

T 400= 400. *1.8 

M= N - 100*M0DE 

BT= PI + P2 

IF (NTAG(BT) .NE.O) GO TO 333 
BWD=BTD / 62.427 

GOTO (1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21,22> t M 

333 BTD= BT 
BTX= BT 

IF (M.NE.6) GO TO 334 
B rCVGL= BT 
BTCO= BT 
RETURN 

334 IF (M.NE.7) GO TO 335 
BTCOND=BT 

335 RETURN 

1 IFIP1.GT.T400) GO TO 101 

BT = 62.427 * BWOENS ( P 1 / 1 . B ^ P2/1^.696) 

GO TO 77 

101 3T = 62.427 * HTDENS ( P 1 / 1 . 8 , P2/ 14. 696 ) 

GO 10 77 

2 IF (P1.GT.T400) GO TO 102 

BT= 1.79882 * B WENHY ( P 1 / 1 . 8 , P2/14.696) 

GO TO 77 

102 BT= 1.79882 * HTENHY (Pl/1.8, P2/14.696) 

GO TO 77 

3 IF (P1.GT.T400) GO TO 103 

BT= .99983 * BWERPY (Pl/1.8, P2/14.696) 

GO TO 77 

103 BT= .99983 * HTERPY (Pl/1.8, P2/14.696) 

GO TO 77 

4 H= 1.79882 * B WENHY (400. , P2/14.696) 

IF (Pl.GT.H) GO TO 104 

BT= 1.8 * BWTOFH (Pl/1. 79882, P2/14.696) 

GO TO 77 

104 BT- 1.8 * HTTOFH ( PI / 1 . 79882, P2/14.696) 

GO TO 77 

5 S= .99983 * BWERPY (400. ,P2/14.696 ) 

IF (Pl.GT.S) GO TO 105 

BT- 1.8 * BWTOFSt PI/. 99983, P2/14.696 ) 

GO TO 77 

105 BT= 1.8 * HTTOFS ( P 1 / . 9998 3 , P2 / 14 . 696 ) 

GO TO 77 

6 IF (Pl.GT.T400) GO TO 106 

BT = .99983 * BWCPGL (Pl/1.8, P2/14.696) 

BTC VGL= BWCVGL*. 99983 
BTCO= BWCO / 2.54/12. 

GO TO 77 

106 BT= .99983 * HTCPGL (Pl/1.8, P2/14.696) 

BTCVGL= BWCVGL*. 99983 

B T C 0 = BWCO / 2.54/12. 

GO TO 77 

7 BT = 1./14.R82 * BWVISC (Pl/1.8, P2/14.696) 

BTCOND* BWCOND # .067163 
GO 10 77 

IF(P1.GT.T400)GO TO 108 
BT = 14.696 * BWPRES(Pl/1.8,P2/62.427) 

GO TO 77 
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108 BT=14.696 * HTPRES(P1/1.8, P2/62.427) 

GO TO 77 

9 P=BWPRES(T 400 »Pl/62*427) 

IF ( P2 / 1 4. 696 • GT • P) GO TO 109 

BT- 1.8 * BWTEMP(P1/62.427,P2/14.696) 

GO TO 77 

109 BT = 1 • 8 * HTTEMP (Pl/62.^27, P2/14.696) 

GO TO 77 

10 14.696 * BWPSAT (Pl/1.8) 

GO TO 77 

11 BT= 1.8 * BWTSAT (Pl/14.696) 

GO TO 77 

12 BT - 62.427 * 2.01572 * ROKSL (Pl/1.8) 

GO TO 77 

13 BT - 62.427 # BWROSV (Pl/1.8) 

GO TO 77 

14 IF (P1.GT.T400) GO TO 114 
BT= .99983 * BWCPIG (Pl/1.8) 

GO TO 77 

114 B T= .99983 * HTCPIG (Pl/1.8) 

GO TO 77 

15 PSAT = BWPSAT (PI / 1.8) + .001 

BT = .99983 * BWCPGL (PI / 1.8 t PSAT) 

GO TO 77 

16 P S A T = BWPSAT (Pl/1.8) -.001 

BT= .99983 * BWCPGL (PI / 1.8, PSAT) 

B TCVGL = HWCVGL*. 99983 
GO TO 77 

17 IF (P1.GT.T400) GO TO 117 

BT = 1.79882 * BWHIDG (Pl/1.8) 

GO TO 77 

117 B T= 1.79882 4 HTHIDG (Pl/1.8) 

GO TO 77 

18 BT= 1.79882 * B WHS TL (PI / 1.8) 

GO 10 77 

19 PSAT = BWPSAT (PI / 1.8) + .001 
B T = 1.79882 * BWHSV (PSAT) 

GO TO 77 

20 IF (P1.GT.T400) GO TO 120 

BT = .99983 * BWSIDG (PI / 1.8) 

GO TO 77 

120 B T = .99983 * HTSIDG (Pl/1.8) 

GO 10 77 

21 BT= .99983 * BWSSTL (Pl/1.8) 

GO TO 77 

22 PSAT = BWPSAT (PI / 1.8) 

BT = .99983 * BWSSV (PSAT) 

77 BTD= BWD * 62.427 
BTX= BWX 
RETURN 
END 

$ I BFTC B I DECK 

B I FUNCTION USED FOR SYSTEM INTERNATIONAL UNITS 

PARA OR NORMAL COMPOSITION ALL PROPERTIES N=1 - 22 PARA 
AND N = 2 0 1 - 222 NORMAL HYDROGEN. PRESSURES NEWTON/SO METER 
AND TEMPERATURES DEGREES KELVIN 
FUNCTION B I ( N * P 1 » P 2 ) 

COMMON/ FIXED/ MODE t ITERS 
COMMON/ BWSONV/BWCO 

COMMON/ B IDXCK/BIDtBI X , B I C VGL , B I COND , B I CO 
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COMMON/BWDXCK/BWD ? BWX t BWCVGL,BWCOND 
DOUBLE PRECISION B ID* B I X , B I C VGL , B I COND* B I CO 
DOUBLE PRECISION BWD ,BWX , BWCVGL ,BWCOND 
DOUBLE PRECISION BWCO 
NODE = N / 100 
T 400 = 400 * 

M = N - 100*MODE 
BI = PI + P2 

IF ( NT AG I B I ) .NE. 0) GO TO 333 
BHD = BID / 1000. 

GO TO < 1,2 ,3,4 1 5,6,7, 8,9 ,10 ,11,12, 13 r 14, 15, 16 ,17,18 1 19,20 ,21,22 ) 

333 BID = BI 
BIX = BI 

IF (M.ME.6) GO TO 334 
B I C VGL = BI 
BICO = BI 
RETURN 

334 IF (M.NE.7) GO TO 335 
BICOND - BI 

335 RETURN 

1 IF (P1.GT.T400) GO TO 101 

BI = 1000. * BWDENS ( P 1 , P2 / 1 01 3 2 5 . I 
GO TO 77 

101 BI = 1000. * HTDENS { P 1 y P2 / 1 0 1 32 5 . ) 

GO 10 77 

2 I E. (P1.GT.T400) GO TO 102 

BI = 4184. * BWENHV (PI, P2/1 01325. ) 

GO TO 77 

102 BI = 4184. * H1ENHY(P1,P2/101325.) 

GO TO 77 

3 IF (P1.GT.T400) GO TO 103 

BI = 4184. * BWERPY ( P 1 , P 2 / 1 0 1 32 5 . ) 

GO TO 77 

103 BI = 4184. * HTERPY ( P 1 , P2 / 1 0 132 5 . ) 

GO TO 77 

4 H = 4184. * BWFNHY(400.,P2/101325.) 

IF (Pl.GT.H) GO TO 104 

BI = BWTOFH t P 1 / 4 1 84 • , P2/101325.) 

GQ TO 77 

104 BI = H1TUFH(P1/4184.,P2/101325.) 

GO TO 77 

5 S = 4184. * BWERPY ( 400 . , P2 / 1 0 1 3? 5 . ) 

IF (Pl.GT.S) GO TO 105 

BI = B.'-'lOFSt Pl/4184. , P 2 / 101325. ) 

GO TO 77 

105 BI = FITTOFSI Pl/4184* *P2 / 101325. ) 

GO TO 77' 

6 IF (PI. GT. 1400) GO TO 106 

BI = 4184. * BWCPGL ( PI , P2/ 101325 . ) 

BICVGL = B WC VGL * 4184. 

BICO = BWCO / 100. 

GO TO 77 

306 BI = 4184. * HTCPGL ( P 1 , P2/ 10 1325. ) 

BICVGL = BWCVGL * 4184. 

BICO = BWCO / 100. 

GO TO 77 

7 BI = .1 * BWVISC ( PI ,P2/1 01325. ) 

BICOND = BWCOND * 418.4 

GO TO 77 

8 IF tPl.GT.T400) GO TO 108 
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in = 101325. * BWPRES(P1,P2/1000. ) 
GO TO 77 

108 B I = 101325.* HTPRES(Pl y P2/1000. ) 
GO TO 77 

9 P = 8WPRES ( T400f PI/ 1000 • > 

IF ( P2/101325..GT.P) GO TO 109 
B I = RWTEMP ( Pl/1000. ,P2/101325. ) 
GO TO 77 

109 B I = HTTEMP(Pl/1000. f P2/101325.) 

GO TO 77 

10 B I = 101325. * BWPSAT { PI ) 

GO TO 77 

11 B I = BWTSAT (Pl/101325.) 

GO TO 77 

12 B I = 1000. * 2.01572 * ROKSL(Pl) 

GO TO 77 

13 B I = 1000. * BWROSV (PI) 

GO TO 77 

14 IF (P1.GT.T400) GO TO 114 
B I = 4184.*BWCPIG (PI) 

GO TO 77 

114 B I = 4184. * HTCPIG(P1 ) 

GO TO 77 

15 PSAT = B WPS AT (PI) + .001 

B I = 4184. * BWCPGL ( PI tPSAT ) 

B I C V G L = 4184.* BWCVGL 
Gf) TO 77 

16 PSAT = BWPSAT (PD- .001 

B I = 4184. * BWCPGL(P1,PSAT> 

BICVGL = 4 1 84. *B WC VGL 
GO TO 77 

17 IF (P1.GT.T400) GO TO 117 
B I = 4184.* RWH I DG ( PI ) 

GO TO 77 

117 B I = 4184. * HTHIDG(Pl) 

GO TO 77 

18 B I = 4184. * BWHSTL(Pl) 

GO TO 77 

19 PSAT = BWPSAT (PI) + .001 
B I = 4184.* BWHSV(PSAT) 

GO TO 77 

20 IF (P1.GT.T400) GO TO 120 
B I = 4184.* BWSIDG1 PI ) 

GO TO 77 

120 B I = 4184. * HT S I DG (PI) 

GO TO 77 

21 B I = 4184.* BWSSTL(Pl) 

GO TO 77 

22 PSAT = BWPSAT (PI) 

B I = 4184.* BWSSV(PSAT) 

77 BID = BWD * 1000. 

BIX = BWX 

RETURN 

END 

$ I BFTC NTAGS DECK 

FUNCTION NTAG(BWS) 

NTAG=0 

RETURN 

END 



APPENDIX C 


DERIVATIONS OF THERMODYNAMIC AND TRANSPORT EQUATIONS 

Equations 

In calculating thermodynamic functions such as enthalpy H and entropy S, or spe- 
cific heats C or C , at either constant volume or pressure, it is desirable to avoid 
v p 

tabular interpolation or curve-fitting each function in two dimensions, but rather to relate 
these functions to the equation of state. To do this requires only that the equation of state 
be differentiable and integrable in the domain of interest. Furthermore, if these opera- 
tions can be performed in closed form, the amount of necessary calculation may be 
greatly reduced. The following derivations demonstrate how to relate H and S to the 
equation of state. 

In general, any state parameter is determined from the specification of any two 
other parameters. Therefore, if S is expressed as S(p,T), entropy can then be for- 
mally written as 


dS = 



dp 


From Maxwell's equation, 




V 


and 


Hence, 


= 1- /9S\ = _-l /as\ _ -1 /9P\ 

" v' W T " p 2 W T " p 2 W, 





dp 


From the second law of thermodynamics, dH = T dS + (dP/p). Setting dP/p equal 
to d(P/p) + P/p^ dp yields 
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dH = T(— 1 dT+df— ] + 

VaT 


_p _ jr (M 

p 2 "p 2 ^ t 'pj 


dp 


Note that along an isotherm the difference between the enthalpy or entropy, H and S, 
at density p, and the values of H Q and S Q at an arbitrary p Q is 


(h p - V T 


_ P(T, p) _ P < T >Po> + 


/ 

*/n 


P _ T_ fdP 

- 2 'p 2 w, p 


dp 


and 


< S P - S P0> = 


n 


9P1 


PJ 


dp 


To avoid the difficulties inherent in integrating through the vapor dome, the problem 
can be subdivided into finding appropriate values of p Q and corresponding expressions 
for H T and S p T in the vapor region and in the compressed liquid region. 

In the vapor region, the problem may be simplified further by calculating 

(H - ' H p,id } and < S n - S n id)_- 


P, id 7 


Since P is equal to pRT for an ideal gas, 


(H . H - H n , H ) = RT - RT + 

P, id P, O, id rp 


/ 


P _ JP 

2 2 y 
\P P /x 


dp = 0 


PfT o) ^C^jPq) 

(H p " H p, id^ T " < H p, o " H p, o, id^ T = ~~ p 


/ 


p t /ap 


J> 2 P 2W PJ T 


dp 


Since it is desirable, as the limit p Q - 0, that the properties of the real gas merge 
with those of the ideal gas and since calculating procedures demand that this be true, it 
is justifiable to use the following relation: 
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lim 

Per 0 


(H-H id ) =^p)_ RT+ f P..JT^SP 


P,T p 


- r 1 

Jo \p 


2 p 2 ^ P 


dp 


If Hj d is then expressed as a function of temperature only and, in the usual manner, 
set equal to H- d (T Q ) + J' Cp(T)dT, the following expression is obtained for H in the 


is 

vapor region: 


H(p, T) = H. 


id< T o> + f 

c/r 


C (T)dT + — - RT + / 

P P / I 2 2 \3T/ 

'T Lp p p 

o 


dp 


The ideal gas entropy deviation at constant temperature is 


( S p, id " S p, o, id^ T 


✓>P> id 
»/> . o. id 1 


-1 /3P1 


L4,o, id p* w p 


dp = 


/ P, id . 

Ji 

, id 


^5) dT 


or 


( S p S p,id^ T " ^ S p,o S p,o, id^ T 


3 


_0P | 


dp + 


*/p, o, id 


/ p,id 
. id 




In the limit as density approaches zero, one can arbitrarily set 


lim(S„ „ - S . ,) = 0 
v P, o p, o, id' 


Thus, 


< s - s id> 


P,T 


3 


_0Pj 

^9T/ 


dp + 


/ p, id 


5 «*• 


Since 
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/ p, id 


rR) 


dp = 




s ^ 


and 


/ 

"P, id 


[— ) dp = R In -P— 
\P/ T P ic j 


where p- d = P/RT. Hence, 


f f-) dp = R In P5! = R In - 
/ \P/ r p P Z 

4, id T 


Therefore, 


(S - S id ) = f [5 - -L (™) 

P ’ T / L p p 2 W 


dp - R In — 
Z 


If S^ d is then 


T 

expressed as equal to S id (T o ) + f C (T)/T(T) dT - R In P where 

y rp 


S id (T Q ) is the ideal gas entr °Py at 1 atmosphere (1. 0132 5xl0 b N/m 2 ); S T in the vapor 
region is given by the expression 


S T,p ~ S id^ T 


t/rc\ 


(T) p 

dT - R In — + 
T Z 


Jo L p p 2 w 


dp 


In the compressed liquid region, the calculation of enthalpy and entropy can be simplified 
by letting p Q be equal to the saturated liquid density for the given temperature. With 
that condition, the enthalpy and entropy differences may be expressed as 
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H(p,T) - H sZ (T) = 


P P sat 


(T) 


P P 


si 


f 

J Psl 


_P _ T_ fdp\ 

2 2 \dT/ 

P P N /r 


dp 


ujt 


and 


S(p,T) - S gZ (T) = 


r fj. /ap\ 

/ 2 \3T/ 

JPei Ip 'pa 


dp 


To calculate Hp gat and Sp gat j. , use can be made of the following properties: 


H SZ (T) - H_(T) = P flat CT) 


sv Vi/ " sat v 


ry 

5 s l P sv/ J | p 
H sv 


p _ /ap\ 

2 ' p 2 W/ p 


dp 


Since the saturation pressure is not a function of density variation in the interval under 
the vapor dome, 


Therefore, 


dP-/^ dT + ^ dp = (^ dT 
\3T/p \3P/ T \9T/ p 



Hence, 


For entropy, 


V T > - H sv< T > = T 
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From these arguments, it is evident that 


- s sv < T > 



Since both H gv (T) and S gv (T) can be calculated as deviations from ideal gas conditions, 
the calculation of H g ^(T) and S g ^(T) would be relatively straightforward. 

The calculations of H g ^(T) and S g ^(T) can be simplified further by considering the 
specific heat of the saturated liquid C g ^(T). For parahydrogen this specific heat has 
been calculated and can be expressed as a function of temperature. Since this function is 
integrable in closed form, the change in entropy along the saturated liquid line may be ex- 
pressed as 





c (T) 

— ~ dT + s W< T o> 


where S_.(T ) is the entropy of the saturated liquid at temperature T . Since from the 

06 O O 

second law of thermodynamics, 


dH = T dS + — 


P 


H 


si 



C s i(T)*T + H gZ (T o ) 



dP 

P 


where H g ^(T Q ) is the enthalpy of the saturated liquid at 
calculated from the saturation boundary equations. 


T and 
o 




dP/p 


may be 



Transport Properties 


The calculation of real gas values of viscosity and thermal conductivity coefficients 
is complicated by the incompleteness of both theory and data. For ideal gas conditions, 
these coefficients are calculated from either theoretical relations (refs. 9 and 10) or em- 
pirical correlations (refs. 5 and 7). 

Attempts have been made to apply real gas corrections by theoretical approaches 
such as that of Enskog (refs. 5 and 11). This approach relates functions of the state 
equation to concepts such as thermal motion and intermolecular forces necessary to de- 
scribe transport properties. The method has met with little success, except at relatively 
low densities. 

Graphical correlations (refs. 9 and 10) have also been developed to represent an 
’’excess function. ” This function is defined as the isothermal increase in viscosity or 
thermal conductivity from the low density limit to the real gas density. These correla- 
tions are based on the hypothesis that the excess function should only be a function of den- 
sity. However, recent studies (ref. 6) have shown that, at least for viscosity, this as- 
sumption is not valid for high densities (i. e. , densities twice critical). A correlation for 
excess viscosity for both temperature and density effects has been developed by Diller 
(ref. 6). The correlation is used to calculate real gas viscosities in this report. It is 
conjectured that further studies will show a similar temperature dependence for thermal 
conductivity. However, in the absence of available data, the excess thermal conductivity 
is at present calculated only as a function of density. 

In all cases, the real gas transport properties are calculated as an ideal gas value 
plus an excess function, as shown in the following equations: 


M - M 0 ( T) + A(p)exp 


B(p) 


Mom= w 3/ !^y 

(T + a 2v )(T + a, v ) 


3v J 


A (p) = exp[b lv + In (p) + b 2v P 3 ^ 2 + b 3y exp(b 4y p)J 


B CP) = b 5v + ^v^ + b 7vP 3/2 + b 8v ^Pfrgvh 3 ) 


K - K q (T) + KjCp) 
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i 



a Ok + a lk T + C p, id (T ^ a 2k + a 3k T) 


I 

K q ( t ) = 



+ a 4k\ 
T / 


K l(p) = b Qk + b lk p + b 2k p 2 + b 2k p 3 + b 4k p 4 + b 5k p 5 + b gk p 6 + b 7k p 7 + b 8k p 8 
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Figure 1. - Region-finding diagrams for pressure-temperature-density input combinations. 





Temperature, T, °K 

Figure 3. - Specific heat of parahydrogen computer plot from library of functions 
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Figure 4. - Density of para hydrogen (computer plot), 
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Figure 7. - Enthalpy of parahydrogen (computer plot). 






1.0x10 s 



1400 


30 40 50 60 70 

Temperature, T, °K 

Figure 8. - Entropy of parahydrogen {computer plot). 1 * 


62 


NASA-Langley, 1968 8 E-4008 








"The aeronautical and space activities of the United States shall be 
conducted so as to contribute ... to the expansion of human knowl- 
edge of phenomena in the atmosphere and space. The Administration 
shall provide for the widest practicable and appropriate dissemination 
of information concerning its activities and the results thereof .” 

— National Aeronautics and Space Act of 1958 


NASA SCIENTIFIC AND TECHNICAL PUBLICATIONS 


TECHNICAL REPORTS: Scientific and technical information considered 
important, complete, and a lasting contribution to existing knowledge. 

TECHNICAL NOTES: Information less broad in scope but nevertheless of 
importance as a contribution to existing knowledge. 

TECHNICAL MEMORANDUMS: Information receiving limited distribu- 
tion because of preliminary data, security classification, or other reasons. 

CONTRACTOR REPORTS: Scientific and technical information generated 

under a NASA contract or grant and considered an important contribution to 
existing knowledge. 

TECHNICAL TRANSLATIONS: Information published in a foreign 
language considered to merit NASA distribution in English. 

SPECIAL PUBLICATIONS: Information derived from or of value to NASA 

activities. Publications include conference proceedings, monographs, data 
compilations, handbooks, sourcebooks, and special bibliographies. 

TECHNOLOGY UTILIZATION PUBLICATIONS: Information on tech- 
nology used by NASA that may be of particular interest in commercial and other 
non-aerospace applications. Publications include Tech Briefs, Technology 
Utilization Reports and Notes, and Technology Surveys. 


Details on the availability of these publications may be obtained from: 


SCIENTIFIC AND TECHNICAL INFORMATION DIVISION 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 

Washington/ D.C. 20546 



